Effects of disorder on electron transport in 
arrays of quantum dots 
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Using analytical and numerical methods, we investigate the zero-temperature transport of elec- 
trons in a model of quantum dot arrays with a disordered background potential, where the electrons 
incoherently tunnel between the dots. One effect of the disorder is that conduction through the 
array is possible only for voltages across the array that exceed a critical voltage Vt- We investigate 
the behavior of arrays in three voltage regimes: below the critical voltage, above, but arbitrarily 
close to, the critical voltage, and further above the critical voltage. For voltages less than Vr, we find 
that the features of the invasion of charge onto the array depend on whether the dots have uniform 
or varying capacitances. We compute the first conduction path at voltages just above Vt using a 
transfer-matrix style algorithm. Though only the first path can be studied using this technique, it 
can be used to elucidate the important energy and length scales. We find that the geometrical struc- 
ture of the first conducting path is essentially unaffected by the addition of capacitive or tunneling 
resistance disorder. We also investigate the effects of this added disorder to transport further above 
the threshold. We find that qualitative behavior is dominated by the presence of the background 
potential, rather than capacitive or tunneling disorder, at least as long as these additional disorders 
do not have an extremely broad distribution. We use finite size scaling analysis to explore the non- 
linear current-voltage relationship near Vr- The scaling of the current / near Vr, I ^ (V — Vr)^, 
gives similar values for the effective exponent /? for all varieties of tunneling and capacitive disorder, 
when the current is computed for voltages within a few percent of threshold. We do note that the 
value of /3 near the transition is not converged at this distance from threshold and difficulties in 
obtaining its value in the V \Vt limit. 



I. INTRODUCTION 



It is now possible to engineer arrays of nanoparticles (1- 
2 nm in diameter) in various geometrical configurations^ 
and to lithographically fabricate arrays of low capaci- 
tance islands separated by tunnel junctions with reason- 
able control over array parameters-^. Such quantum dot 
arrays (QDA) have been the subject of intense investiga- 
tion recentlySii. In spite of the relatively well controlled 
properties of these arrays however, there are limitations 
on the homogeneity of such systems. Disorder at the 
sub-micron length scales arises due to a variety of rea- 
sons, is inevitable and significantly influences the prop- 
erties of these otherwise well ordered arrays. For ligand 
coated nanoparticles the variation in coating properties 
and separation result in different resistances to electron 
tunneling. As a consequence of the poly-dispersion in 
the sizes of metallic nanoparticles, the charging energies 
of the individual islands differ. Similarly for lithograph- 
ically fabricated tunnel junctions, islands with variable 
charging energies arise due to a dispersion of island sizes 
or due to fluctuating capacitative coupling between dots 
and the underlying gate. Given the pervasiveness of ran- 
dom background charges, nonuniform charging energy 
and fluctuations in tunneling resistance across the array, 
an important focus of the current work is to study the ef- 
fect of these on transport properties of electrons. Due to 
limitations of the fabrication process, it is difficult to con- 
trol the different types of disorder independently, whereas 
this can can be relatively easily addressed by computer 
simulations. 



There are many dynamical systems in which strongly 
interacting particles exhibit collective transport in a ran- 
dom environment^. Even though the underlying micro- 
scopic details are different, systems like the vortex glass 
in type-II superconductors and charge density waves&, 
share some general features in the long wavelength limit, 
e.g., they are both characterized by the presence of a 
well defined threshold force below which the system is 
essentially static and above which the system has a non- 
linear response. Electron transport in disordered QDA 
also provide a useful system to study problems of quali- 
tative similarity. An advantage of QDA is that the pri- 
mary interactions and the fundamental physics are rela- 
tively better understood and arguably under greater ex- 
perimental control. 

Using a combination of analytic and numerical tech- 
niques, we investigate electron transport at zero temper- 
ature in arrays of disordered small capacitance islands 
which are capacitatively uncoupled to their neighbors. 
We use this system both as a model for collective trans- 
port of discrete charges in a random environment and for 
better understanding the role of disorder. By studying 
similar systems, but for different values of parameters, 
different regimes of the collective transport problem can 
be addressed. These regimes may be characterized by 
the relative strengths of disorder, tunneling rates and 
electron-electron interaction. These regimes are accessi- 
ble experimentally too, as arrays can be fabricated with 
varying degree of tunability of the coupling between the 
array elements^i. For example, the model in Refs. 
Q is similar to the model we study in this paper - in 
that offset charge disorder is included, although Gaus- 
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sian distributed as opposed to uniformly distributed - 
but the screening length is assumed infinite. (Transport 
in this regime is believed to belong to the same univer- 
sality class of two dimensional magnetic vortex model in 
disordered superconducting films.) Ref. [loj which inves- 
tigates 2DEG at semiconductor heterointerfaces designed 
to keep the self-capacitance and disorder low - the I-V 
characteristic is better explained by charge soliton injec- 
tion as no threshold voltage is observed. Such systems 
can be used to study two dimensional Coulomb gases 
which undergo charge Kosterlitz-Thouless (KT) transi- 
tions. 



A. Outline 

There have been many experimental and simulation 
papers ^LiSiiiiiii^ (to cite just a few) which have used 
the theory developed in Ref. by Middleton and 
Wingreen (MW). This paper expands on the MW dis- 
cussion of electron transport in disordered arrays. The 
original model and its extension to include other forms 
of disorder are described in the remainder of this section. 
One-dimensional arrays - both below and above thresh- 
old - are discussed in detail in sectional In section UTTI 
results for 2D arrays below the threshold voltage are pre- 
sented including several results not discussed previously. 
As the original MW paper sketched only briefly the con- 
nection between the independent conducting paths and 
the properties of a directed polymer in random media 
(DPRM), a major aim of section llVl and llTTl is to establish 
the connection on a more rigorous basis. In section llVl we 
discuss the morphology and current carrying properties 
of the first conducting path at threshold for QDA. It also 
provides some of the details required to understand the 
non-linear scaling of current (I) with voltage (V), which 
is addressed in section Ivl 

There have been several papers that have used nu- 
merical approaches to investigate transport in arrays 
(both ID and 2D) in the presence of random background 
charges as well as other types of disorde r^^i^^i^^i^^ . Some 
approaches, have used discrete event simulation tech- 
niques to model the individual tunneling events, whereas 
some have explicitly used computed transition rates in a 
master-equation approach. The common aim is to com- 
pute the general TV characteristics, which as a conse- 
quence of the collective behavior of electron tunneling 
is non-trivially dependent on the individual rates. We 
focus on a statistical physics approach to the problem, 
thus laying a theoretical basis for the scaling exponents 
observed experimentally and numerically. 

B. The Model 

The three main energy scales of QD are the charging 
energy (-Ecs), the electron in-a-box energy levels (A) and 
the thermal energy (fcT). As a consequence of the small 



size of these islands and tunnel junctions the capacitance 
involved are in the femto to atto-Farad range, thus the 
charging energy - which is the increase in energy due to 
the addition of a single electron is given by e^/Cs - of 
these islands is large. A characteristic feature of QD is 
the clear separation of internal energy scales A and E^^ . 
An external energy scale (fcT) is set by the temperature 
of interest, which determines the levels that are resolved 
and participate in transport. When Ec^^ kT the role 
of thermal fluctuations can be ignored. Depending upon 
the temperatures of interest, A maybe comparable to 
kT or different; for kT ^ A, the discrete energy level 
spectrum of the QD do not play a role during transport. 
Metallic dots are different from semi-conducting dots by 
the fact that typically the level spacings for metallic dots 
are much smaller compared to other energies. At suffi- 
ciently low temperatures, the scale of which is set by, Ecj, 
^ kT, the addition of a single extra electron to an dot 
increases the dot energy; in spite of the increased energy, 
the dot is stable to thermal energy fluctuations, which in 
turn makes it unfavorable for more electrons to tunnel 
onto the same dot, resulting in its blocking other elec- 
trons onto the dot. This is called the Coulomb blockade 
regime. 

The parameters required to characterize QDA can vary 
over a large range of values and consequently so do the 
properties of QDA. Thus, it is instructive to understand 
the parameter space of QDA in order to appreciate the 
details of the model. The main parameters used to char- 
acterize QDA, as opposed to individual quantum dots 
(QD) are: the tunneling resistance (Rt) which to a first 
approximation is a measure of how well confined the elec- 
trons are on a dot, the inter-dot capacitance (C/) and the 
dot capacitance (Ce) which is a function of the junction, 
gate and self-capacitances. The relative values of Ce 
and C'l are important as it determines the extent of elec- 
trostatic coupling between dots in the array. The exact 
value of Cs depends on the system under consideration. 
For example, typically the self-capacitance of lithograph- 
ically prepared arrays is negligible compared to the other 
capacitances, thus Cs is a function of the dot-gate and 
tunnel junction capacitances (for example in Ref. |l3j | Ce 
= Gg + 4C). For nanoparticles with diameters of a few 
nm, the self-capacitance becomes important and should 
possibly be considered in the computation of C^r^- Ei- 
ther way, Ce still sets the scale for the charging energy. 
Independent of the actual experimental setup considered, 
as long as Cs ^ C/, the dots are considered to be ca- 
pacitatively uncoupled to each other and the electrostatic 
energy is determined by on site interactions only. How- 
ever if Ce is comparable or less than Cj the dots are 
capacitatively coupled. A screening length (A) can be 
thought of the distance (in units of dots) upto which the 
charge on a dot can be felt electrostatically, i.e., distance 
that an excess charge placed on a dot will effect neigh- 
boring dots by polarization. The polarization decreases 
exponentially with A, which in turn decreases with the 
ratio of S^; this is consistent with understanding that 
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there is a stronger screening of the electrons on a dot 
from electrons on adjacent dots as the capacitative cou- 
pling between a dot and the back gate increases. For Cs 

The main modification to the original MW model is the 
introduction of nonuniform dot capacitance Cs and tun- 
neling resistances Rt- The effects of underlying charge 
impurities trapped at the interfaces and in the substrate 
are captured in a random background charge on each 
dot. The effect of background charges is modeled as off- 
set charges on each dot (qi). The offset charge at any 
site is considered to be [0,1 [, as any value outside this 
range will be compensated by electron hopping. Arrays 
with only offset charge disorder are referred to as UC 
(uniform capacitance) systems. The area and capacita- 
tive coupling of an island to the underlying electron gas 
varies from dot-to-dot in an array. These fluctuations in 
the dot-gate and self-capacitance of dots, along with stray 
capacitances are incorporated by assuming a varying dot 
capacitance Cs. As Cs controls the charging energy of 
the dot, a non- uniform Cs results in different charging 
energies Ecj^oi dots. Arrays with both offset charge dis- 
order and a varying Cs, are referred to as DC (disordered 
capacitance) systems. Fluctuations in the tunneling re- 
sistance - either due to varying distance between metal- 
lic dots or the varying material properties of the tunnel 
junction separating the metallic islands between dots - is 
captured by assuming a log-normal distribution of tun- 
neling resistances. Arrays that incorporate variation in 
tunneling resistance as well as offset charge disorder, but 
with a fixed value of Cs , are referred to as RT (resistance 
disorder) systems. 

We assume small metallic islands are separated from 
each other by tunnel junctions of resistance Rt but ca- 
pacitatively coupled to neighboring dots (C/). We as- 
sume a constant capacitance C/ between neighboring 
dots and between the left and right leads and dots adja- 
cent to them. The dots are assumed to be separated from 
an underlying back gate by an insulating layer. Each 
dot is capacitatively coupled to the back gate with a 
capacitance Cs. The leads and back gate are assumed 
to have infinite self capacitance. As a consequence of 
the proximity of the back gate to the dots Cs ^ C/, 
the screening-length is taken to be less than one lattice 
spacing. Consequently the capacitative coupling between 
dots is neglected. 

We will consider arrays where the single-energy levels 
of the dots are essentially a continuum at the Fermi level 
in the strongly Coulomb blockaded {Ec^'^ kT) regime. 
Thus tunneling is between levels determined by E^^ . 
Where a spread in values of Rt is considered, we as- 
sume tunneling resistance between any two dots is still 
sufficiently large to consider electrons localized on a site 
{Rt ^ h/e^). This is the regime of the "orthodox the- 
ory" of single electron tunneling and is applicable for 
both the micron sized lithographically defined SET (e.g., 
metal islands embedded in a substratei^ and separated 
by tunnel junctions or semi-conductor islands separated 



by barriers-'') as well as the 3D metallic grains. According 
to the "orthodox theory" of a tunneling event across a 
tunnel junction, tunneling rates (transition probability 
per unit time) associated with an event are given by. 



where Ai? is the difference in the free energy of the sys- 
tem before and after the tunneling event, Rt is the tun- 
neling resistance of the junction involved in the tunneling 
event, T the temperature and k is the usual Boltzmann 
constant. The kinetic energy gained by the tunneling 
electron is assumed to be dissipated. The value of Rt is 
assumed to be much greater than h/e^. This essentially 
implies that the wavefunction of electrons are localized 
to a single dot which permits the number of electrons on 
any single dot to be treated as a classical variable, ft 
should be pointed out that the orthodox theory is still 
valid for arrays in the limit Ci ^ Cg , but not for dots in 
the other limits of Rt <C Rq and Ec^<^ kT. 

In this limit the energy is all electrostatic and is de- 
termined by a capacitance matrix dj and is represented 
as: 

E = VlQl + VrQr + \Y1 + <l^)C'^Qi + (2) 

where Ql [Qr] are the charges on the left (right) leads, 
which are at voltages Vl {Vr) and C"*-' is the inverse of 
the matrix of capacitances between dots i and j. The 
diagonal elements of Cij are the sum of all capacitances 
associated with a dot and the off-diagonal elements are 
the negative of the inter dot capacitances. Hence for 
a NxN array in the limit of ^ ^ 0, the capacitance 
matrix is a NxN diagonal matrix. 

In the limit of small screening length (less than 1 dot 
spacing) and the presence of offset charge disorder the 
voltage on dot i is given by Vi is {Qi + qi)/Cs- 

At zero temperatures the expression ^ for tunneling 
rates reduces to 

A F 

r = ^e{AE) (3) 

hence a charge may tunnel from dot i to j, only if such 
an event lowers the free energy of the array i.e. 

V, > Vj + e/Cs (4) 

II. ID ARRAYS 

Before attempting to understand the detailed proper- 
ties of two dimensional arrays, we begin by an attempt 
to understand the relatively simpler case of a linear chain 
of quantum dots, as they facilitate an understanding of 
some of the ideas required later. There have been sev- 
eral experiments aimed at understanding the conduc- 
tion properties of essentially one dimensional arrays of 
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nanoparticleaSSiS^. The ability of metallic nanoparticles 
to be patterned using polymers templates makes them 
attractive candidates for potential future self-assembling 
electronic devices. 



A. Uniform Ce: Insulating State 

We start by exploring the tunneling of electrons onto 
the array from the emitter lead. In the zero temperature 
limit, as the capacitance of the leads is assumed to be in- 
finite, electrons can flow onto the array when the voltage 
of the emitter lead (Vl) is equal to or greater than the 
voltage of the leftmost dot as given by Eqn. 0] At this 
applied voltage, an electron cannot tunnel from the left- 
most dot to the next dot, say represented by the index i, 
if dot i has an offset charge impurity qi greater than the 
offset charge impurity q of the leftmost dot. Electrons 
tunnel onto the array only if it is possible to do so with- 
out an increase in the free energy of the system . Thi s is 

In 



no longer possible for the configuration in Fig. l(a 



this configuration the electron residing on the leftmost 
dot is considered to be pinned. 



As shown in Figs. 1(b) and 1(c) the emitter lead volt- 
age has to be increased by at least one unit in order that 
the electrons can overcome the barrier. As the value of 
the emitter lead voltage is successively increased, there 
will be a cascade of electrons tunneling onto the array 
from the lead, until they get progressively pinned and it 



is no longer energetically possible for electrons to pene- 
trate further into the array. The flow of charges onto the 
array at a given emitter lead voltage until they are all 
pinned due to the disorder and thus no further electrons 
can tunnel onto the array constitutes an avalanche. 

There exists a unique value of the emitter lead voltage 
- which depends upon the underlying disorder profile - at 
which electrons wil l be a ble to reach the collector lead for 
the first time (Fig. 1(c) I. This well-defined voltage value 
is referred to as the threshold voltage {Vr). Vr separates 
the conducting phase from an insulating phase. Typi- 
cally in order to reach the collector lead end of an array 
L dots long, an electron will have to overcome ^ upward 
steps. These steps can be can be understood as the av- 
erage number of steps a random walk in ID makes in a 
given direction, thus the mean threshold voltage should 
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FIG. 1: A schematic illustrating the build of charges in a 
ID array as the emitter voltage is progressively increased to 
threshold. 



FIG. 2; In Fig |2(a)| a plot of the scaling of Vr and root mean 
square fluctuations of Vt w ith sy stem size for a ID array when 
the Ce is uniform. In Fig |2(b)| scaling of Vt and root mean 
square fluctuations in Vt with system size for a ID array but 
with non- uniform distribution of Cs- 
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be given by 



(5) 



where the over-bar represents an averaging over disorder 
reahzations. Sample-to-sample fluctuations in the Vt can 
be thought of as the root-mcan-squarc fluctuations of a 
random walk in ID which scales as N^^^, where N is the 
number of steps of the random walk. Hence fluctuations 
in Vt should scale with system size as 



(6) 



The scaling o f both Vt and (t{Vt) with system size, as 
shown in Fig. 2(a) are consistent with the above expla- 
nation. 



B. Uniform Ce: Conducting State 

The threshold voltage represents the lowest voltage at 
which electrons can tunnel across the array, hence for 
Vl > Vt, current flows through the array. For a given 
disorder realization, Vt depends on the number of up- 
steps encountered due to the offset charge impurities. 

If Vl is marginally greater than Vt, so that = 
{Vl — Vt)/Vt <C 1, then the discreteness of charge and 
offset charge impurities play a crucial role in determining 
the current. At these voltages the current is determined 
by the slowest tunneling rate {Tgiow) between any two 
neighboring dots in the array (analogous to the net flow 
of traffic being determined by the bottleneck in the path 
of flow), which on the average is given by where 
L is the number of dots in the ID array. 

This can be understood as follows: (Vl — Vr) represents 
the voltage increment over Vt- In principle the voltage 
drop can be anywhere between between and Vl — Vt 
for a given pair of dots. For an array with L dots there 
are L-l-1 (~ L when L is large) voltage drops (tunneling 
rates), hence the minimum voltage drop across any two 
dots is on the average , which results in a tunneling 

rate given by Eqn. ^ to be Using Vt = we 

get Tsiow = 2RCsVt • ^ ^ e^siow we have 



(7) 



As can be seen from simulation results in Fig. |3(d)' 



the limit of low i> and for large system sizes, the local 
value of the exponent is consistent with 1. Note that the 
for smaller system sizes (L < 1000) the effective exponent 
is quite far away from 1. The fact that chains at least 
larger than 1000 are required is an important observation. 
We will revisit its implication later in the chapter. 

In the opposite regime of a high applied voltage, ^ 
1, the current is determined by the average tunneling 
rate across a pair of dots. This is given by the ^ of 
the average voltage drop across a pair of dots, ^-^7^^ , 



i.e., r = ^ej^j/^ , which gives the same scaling expression 
for the current with v as Eqn. (Q. For values of z/ ~ 
1, a crossover from slow point dominated current linear 
scaling to high applied voltage linear scaling is observed. 



C. Non-Uniform Cs: Insulating State 

As mentioned, the introduction of dots with non- 
uniform Cs leads to an array with dots of different charg- 
ing voltages. For our simulations, we assume Cs to be 
uniformly distributed between 1.0 and a maximum fluc- 
tuation of 2.0. If we attribute the variation in Cs by 
a factor of 2 to fluctuations in the size of the dots, it 
corresponds to a change in a variation in the linear di- 
mension of dots by a factor of V2. The determination of 
the Vt gets complicated by the presence of both offset 
charge disorder and varying charging energies. This is 
illustrated in Fig. |4(b)| where the spacings between volt- 
ages are different for different dots; this is in contrast to 

The incre- 



dots with uniform Cs as shown in Fig. 4(a) 



ment in voltage required in order to tunnel between a 
pair of dots is due to two independent random variables 
- ^ and Pi, where l/Cs is the charging energy of dot i 
and f3i is between and 1 which represents the required 
increment due to the offset charges. Vt can be written as 
the Y,i{[3i / Cy:) , where the summation runs over the num- 
ber of dots in the array, L. Vt can therefore be written 
as L.{l3i) .{-^) , where () represents the average values. 
Hence 



Vt - (l/Cs)i/2, 



(8) 



where ^ is log(2.0) for the assumed maximum value of 
2.0. Although Vt scales as L, similar to the UC arrays, 
cr(VT) behavior is more complicated. An exact analytical 
expression for cr(VT) - which can be derived using the 
expectation values of 1/Cs and 1/Cs^ gives. 



<j{Vt) 



f-^(log(2))^. 



(9) 



Thus, up to leading order (t(Vt) scales as L^/"^. This is 



consistent with our results as can be seen in Fig. 2(b) 

An often used technique to explore the disorder en- 
ergy scale is to study the response of the system on 
changing the boundary condition^^. For disordered ID 
QDA, we change the boundary condition at the right 
lead and study the change in threshold voltage. We de- 
fine AVt(^Vr) as the difference in Vt on changing the 
value of Vr by SVr- Recall that for UC arrays Vt was 
completely determinable by the number of up-steps in 
the offset charge impurities. For the uniform capacitance 
case the response is trivial: a shift in the by SVr - 
where 5Vr is ^ or a multiple thereof - changes Vt by 
the same amount. This is a consequence of the response 
being periodic in voltage. For the ID QDA with disor- 
dered capacitances, however, a shift in Vr by SVr guar- 
antees a change in Vt by SVn only on the average, due 
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FIG. 4: Fig. |4(a)] is a schematic of dot voltages in a ID array of 
dots with uniform Cs- The block heights indicate the increase 
in potential on the addition of an electron. The block heights 
are the same for different dots as Ce is the same. The broken 
line boxes at the top indicate the potential if another electron 
were added. The left arrow indicate the difference in potential 
that will determine the rate of tunneling when the occupation 
number of the leftmost dot is 4 and its right neighbor is 2. 
Similarly the right arrow indicates the rate of tunneling when 
the occupation number of the last dot is 1 and the preceeding 
dot to its left is 2. Fig. |4(b)| is a similar schematic when Ce 
is non-uniform, with Ce for the second, fourth and fifth dots 
different. 
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FIG. 3: I-V curves for ID arr ay wi th uniform capacitances 
drawn on a log-log plot. Fig. |3(b)| is the plot of the local 
slop e of t he current versus u derived from the data in plot 
Fig. |3(a)| The value of the local exponent is computed us- 
ing the values of the current at neighboring values of ly, and 
assigned a value equal to the geometric mean of the two v. 
For clarity small systems have been separated from the largest 
three ID chains simulated. In general at low voltages the cur- 
rent scales linearly with reduced voltage {ly). Flat regions for 
smaller system sizes arise because the corresponding increase 



to a consequence of the threshold voltage being invariant 
to the zero-level of the lead voltages. Specific values of 
AVt depend upon the specific disorder configuration. As 
a consequence of the invariance just mentioned. 



{AVt) = P{AVt + 0)(AyT)Ay^^o + 



(10) 



where P{AVt ^ 0) ( P{AVt = 0) ) is the probability 
that the threshold voltage changes (does not change), 
(AVt) aVt#o the average of the non-zero values of AVt — 
0. Also {AVt)avt=q is 0. 

The response to a change in the right lead voltage can 
be formulated in terms of a ID random walk problem. 
Given the initial and final points of a random walk, one 
can ask what is the probability that a random walk start- 
ing a distance a from the initial point of the original walk. 
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expect that the probabihty of AVr (i.e., survival) 
and the mean of the non-zero AVr ((AVr)AVT#o) should 
scale as L~2 and La respectively. This i s cons istent with 
our numerical results as shown in Fig. 5(a) , although 



there are significant deviations at smaller system sizes. 



D. Non-uniform Ce: Conducting State 

We discussed earlier how for UC arrays in the regime 
of low v, the value of the current is determined by the 
presence of dynamically important slow points. An im- 
portant distinction that arises in DC arrays is that the 
location and value of slow points is less regular. For 
UC arrays, the value of smallest voltage drop - and 
hence the minimal tunneling rate - was bound to in- 
crease every time the emitter lead voltage was incre- 
mented by one unit (-fj-)- Unlike UC arrays, the amount 
by which the emitter fead voltage must be increased in 
order to overcome the slow point does not have a well- 
defined lower bound and varies significantly from sample- 
to-sample and with the value of reduced voltage. As can 
be seen from the presence of plateaus in Fig. El the I-V 
at low v for a single DC array is qualitatively different 
to a sample-averaged I-V. 

At a given value of the mean tunneling rate (F) 
is proportional to the average potential gradient, and is 
given by, 



(11) 



(b) 



FIG. 5: Fig. |5(a)| plots the probability for ID systems with 
non-uniform Cs that that Vt changes when the right lead 
voltage is increased by one unit. The probability decreases 
as L~2. For ID systems with non-uniform Ce, the mean 
value of non-zero AVr scales as the square root of the system 
size. As a consequence of invariance, the mean value of all 
AVt is equal to the value by which the right lead voltage is 
incremented. 



intercepts the original walk before a distance 17 If we as- 
sume that interception with the original walk results in 
annihilation, we can ask of the surviving walks - what is 
the typical separation of the end-point of surviving walks 
from the end-point of the original random walk? It is 
known^^, that the probability of "survival" decreases as 
l^/^ and the typical separation scales as l^^^ (the square 
root of the mean standard deviation of a I step random 
walk). The mapping to the random walk problem is car- 
ried out by considering Vr as the origin of the walk, the 
potential of each dot at threshold (minus the gradient) to 
be the positions of the original random walk and finally 
SVfi as the distance a of the initial point of the second 
random walk from the original random walk. Thus, we 



Based upon the relative values of the F and the typical 
maximum fluctuations from F, we can categorize the ap- 
plied voltage into three regimes. These regimes are: (i) 
when the maximum fluctuation are larger then the mean 
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FIG. 6: I-V characteristic for a single sample ID QDA with 
disordered Ce- Unlike UC arrays, at small reduced voltages 
an the current value remains constant over a range of v val- 
ues - hence the observed plateau (s). This happens when the 
smallest tunneling rate does not change in spite of increasing 
v. When the rate at the slow point does change, however, the 
value of the current jumps resulting in the step like features. 



8 



tunneling rate; (ii) when the maximum fluctuations are 
of the same value as the average gradient, and (iii) when 
the maximum fluctuations are much less then the average 
gradient. 

In regime (iii) the fluctuations about the mean gra- 
dient can be ignored; they no longer influence the cur- 
rent value and consequently the current is given by the 
Eqn. The average potential profile at any given dot 
can be computed using the average potential gradient and 
the distance of the dot from the boundary. By subtract- 
ing the dot potentials from the averaged potential profile 
at each site, we can calculate the fluctuations, and thus 
the roughness of the voltage surface. In regime (i) the 
roughness of the voltage surface scales as L2. Assum- 
ing Gaussian distribution of the fluctuations, the mean 
value of the maximum of N variables from a distribu- 
tion with mean fi and standard deviation a is given by /x 
+ a^/a log N"^^. The deviat ion from the mean tunneling 
rate is maximum at the slow point. Given that the slow 
point can occur anywhere between any two dots (i.e., 
and L), the typical value of the maximum fluctuation 

scales as fi + a log(-|-). The increment in voltage /S.v 

should thus be greater than the maximum barrier (fluc- 
tuation) in order that the slow point be overcome, i.e., 
an additional electron flows over the slow point which in 
turn will result in an increase in current. Therefore the 
probability that a change by Aj/ will overcome a slow 
point is given by the probability that the typical maxi- 
mum fiuctuation is less than Az/. As /S.v ^ L, P(Az^ > 

which approaches 1 
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Li/2^alog(L)' 



typical maximum) 

as L gets larger. 

If step like features persist in the I-V for single sam- 
ples, then given the sample-to-sample fluctuations in the 
location of the plateaus, the sample avera ged I -V curve 
will be more or less flat. As seen in Fig. |7(a)| there is 
a voltage upto which the averaged current is more or 
less static. The value of this voltage decreases with in- 
creasing system size - consistent with the arguments that 
the same increase in v is more likely to result in a slow 
point being overcome as L gets larger. In spite of the 
irregular change in the value of the minimum rate, the 
average value of the minimum voltage drop across any 
two dots remains ii^ii. Hence the average value of 
the minimum rate remains as before - ^''nY^ , which im- 

eHL ' 

plies that once the "static current regime" is overcome 
the current should scale linearly with voltage. Thus in 
spite of the introduction of variable Cs, current scales 
linearly with in regimes (i) and (iii), similar to UC ar- 
rays. The crossover from linear scaling in regime (i) to 
(iii) - corresponds to regime (ii) and is more complicated 
to understand analytically. 



In Fig. 7(a) and Fig. 7(b) we plot the current and lo- 
cal exponent values for the largest systems (L > 500). 



As shown in the plot of effective exponents (Fig. 7(a) I, 
numerically we flnd /3= 0.85 ± 0.02 in the low voltage 
regime (0.05 < ly <0.5) and 1.2 ±0.05 in the high voltage 
regime {ly > 1.0). Similar exponent values are found for 
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FIG. 7: Analogous to Fig. |21 the I-V curves and Piocai for 
ID arrays with disordered Ce . The major difference is that 
the value of the effective exponent, even for the largest ID 
arrays, does not appear to plateau at 1.0, but seem to flatten 
out at 0.85. 



system sizes less than L=500 (not shown). 



III. 2D ARRAYS: INSULATING STATE 

We saw in the previous section, how for one- 
dimensional arrays, charge flowed onto the array from 
the emitter lead till it was energetically favorable. In 
this section we will attempt to develop an understanding 
of the progressive build up of charge in two-dimensional 
arrays, as the emitter lead voltage (Vl) is increased; the 
tunneling of charge is still governed by Eqn. ^ The flow 
of charge onto the array can thus be viewed as lowering 
the energy. Such relaxation of charges so as to lower the 
system energy, is analogous to several different systems 
where the system reaches a lower energy via a series of 
avalanchea^*S&. 
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The threshold voltage is the minimal emitter lead volt- 
age possible such that when electrons tunneling onto the 
array from the emitter lead have sufhcient potential to 
overcome the disorder barriers and reach the collector 
lead. Given a disorder configuration it is not trivial to 
determine the minimal voltage for 2D arrays. A naive 
approach might be to think of the LxW array as W, ID 
arrays of L dots each; trivially compute the "threshold" 
voltage for each of the ID arrays and then find the min- 
imum. The computed minimum would still probably be 
overestimating the true threshold voltage. Determining 
the threshold voltage can be formulated as an optimiza- 
tion problem, but motivated by the aim of understanding 
the physical buildup of charge in QDA, we take a different 
approach. For a given emitter voltage, we add charges 
till a meta-stable insulating state is reached; then the 
emitter lead voltage is progressively increased, building 
up charges until an insulating state no longer exists. The 
value of the emitter lead voltage at which electrons first 
tunnel onto the collector lead is our computed threshold 
voltage. 



A. Avalanches 

As briefly mentioned earlier that an avalanche at a 
given voltage is the flow of charge onto the array un- 
til the flow is arrested by disorder. Avalanches in QDA 
are qualitatively similar to those found in other systems 
with collective elastic transport. Some well studied ex- 
amples are vortex flow in disordered superconductorsSi 
and the avalanches when an interface like a CDW moves 
in quenched disordered system a^^i^^ . For 2D arrays, the 
location where charges tunnel in a given avalanche, helps 
develop the notion of connected elastic domains - basins. 

We define qi{V]^) as the charge of site i before the 
emitter lead voltage is incremented to Vl and qi{V^) as 
the charge of site i after the emitter lead voltage has been 
incremented to Vl ■ The physical size A of an avalanche 
is the number of sites where qi{V]^) ^ qi{Vj^), and the 

volume is J2i Qii^h) ^ Qii^L)- If '^6 set n{A, Vl) to be 
the number of avalanches between size A and A + dA, 
at an emitter lead voltage of Vl and define N{A) — 
J^^ n{A, V)dV , then N{A) can be thought of as the 
number of such avalanches that occur in going from a 
Vl= to = Vt. 

We explore the distribution of avalanche sizes for 
square samples (LxL). The size of an avalanche A can 
also be thought of as the "surface area" - which is equal 
to the number of dots that electrons tunnel onto during 
an avalanche at a given Vl- As the size of avalanches 
vary over several orders of magnitude - starting with 
avalanches of size 1 to system spanning avalanches - we 
use logarithmic bin sizes. Logarithmic binning is natu- 
ral for exploring power laws and reduces fluctuations in 
plots. 

Using standard finite size scaling, we conjecture a scal- 
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FIG. 8: Scaling collapse for the distribution of avalanche sizes 
N(A), for symmetrical systems with uniform Ce is plotted in 
Fig. |8(a)| The exponent a, gives the typical size of the largest 
avalanches as L" ; avalanches of siz es greater than L" become 
increasingly improbable. Fig. |8(b)| shows the collapse of data 
for the mean linear size of avalanches with size between A 
and A+dA for systems with uniform Ce- From the scaling 
collapse we estimate df to have a value of (= p/a) = 1.5. 

ing form for N{A) to be of the type: 

N{A) = L^NiA/L""), (12) 

where N{x) scales as x'^ in the limit of a; ^ 1 and N{x) 
approaches a constant in the limit of a; ~ 1, where L 
is the length of the system. The exponents a and b are 
determined to be those exponents for which a scaling 
plot of N{A)/L'' vs A/L"" yields a single scaling function 
N{A/L°'). The two exponents are not independent and 
can be shown to be related by the relation a -I- 6 = 34Si In 
addition to the two exponents a and 6, a third exponent 
c, can be used to make the curve flat in the regime where 
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defined in the scaling function: 

1{A) = L^LiA/LP), 



(13) 



to have values consistent with |and 1 respectively. We 
get a collapse to l/L"^ ^ {A/L'')'^ for all system sizes L, 
thus I A'^ and the relation constraining the exponents 
is therefore cr = pn ov n = ajp. We know that the 
A^ V^i ^ where c?/ is defined to be the fractal dimension, 
from which we get n, = thus df = p/o. From the 

computed values of p and tr, d/ works out to be 1.5. This 
is consistent with the conclusions from the distribution 
of avalanche sizes. 

Finally, we have also investigated the avalanche struc- 
ture using the radius of gyration {Rg) of avalanches, 
which is defined in the usual way as: 



10 



0.1 



0.01 



512x512 
343x343 
216x216 
125x125 
64x64 
27x27 



-4 



p= 1.0 
c=0.67 



0.001 0.01 



0.1 



10 



100 1000 



(b) 



FIG. 9: The data collapse for distribution of avalanche sizes 
for symmetrical arrays with non-uniform Ct. is plotted in 
Fig. |9(a)| whilst the collapse of data for the mean linear size 
of avalanches with size between A and A-l-dA for sy mmet- 
rical systems with non-uniform Ce is plotted in Fig. |9(b)| 
From the range of exponents for which the scaling collapse is 
acceptable we get the d/ (= p/cr) = 1.5 ± for avalanches. 



A < L° , which for square systems of length and width L 
is related to the other two exponents by 6 — ac = 2 4^ As 
there are two constraints for the three scaling exponents, 
we get only one independent exponent from the scaling 
collapse of th e dist ribution of the sizes of avalanches. As 
shown in Fig. 8(a) the collapse of data to a single scaling 
function is satisfactory, which indicates that the dimen- 
sion of the avalanches is | . The typical size of the largest 
avalanches is given by . To study further the morphol- 



ogy of the avalanches, we compute the the mean of the 
maximum length of avalanches with sizes between A and 
A-|-dA. We collapse the data as shown in Fig. 



8(b) 



onto 



R 2 _ E(?-. 

Kg - ■ 
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N 
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a single curve and determine that exponents a and p, 



and study the scaling of the mean Rg for avalanches of 
sizes between A and A + dA. Numerical evidenca^S indi- 
cates that the scaling of the area with Rg is similar to the 
fractal dimension of the avalanches, which implies that 
the avalanche morphology is compact, i.e., does not have 
any significant holes. 

We have investigated avalanche structure using three 
ways and the results of all three are consistent with the 
hypothesis that typical avalanches are compact with di- 
mension of |. For systems with uniform Cs, the sequence 
of dots at which avalanches originate is periodic in left 
lead voltage (with periodicity ^)- We can thus think 

of "basins" of dimension | evolving as charge flows into 
the array, with some basins growing at the expense of 
others. In general the basin structure is not isotropic, 
as they have a preferred growth direction and the linear 
size in the direction transverse to this preferred direc- 
tion grows only as where I is the linear extent in the 
preferred direction. Thus in a square samples of length 
and width L, there are approximately Nb{l) ^ L/la in- 
dependent regions of activity, where Nbil) is defined as 
the number of basins at a distance I from the left lead. 
Hence to increase the chances of having the large basins 
that scale with L, we simulated systems of length L and 
width a multiple of Ls (width = Nb Ls). 

Similar to square samples, exponents in the scaling col- 
lapse for the distri bution of avalanche sizes for systems 
of size LxLS (Fig. 10(a) I, are not independent but con- 
strained by two relations. 

Given that the total number of dots is the sum of 
the product N{A)A scales as La rather than L^ °; hence 
a + b = |. Also the number of avalanches in the bin [A, 
A-(-dA] scales as , so — ac = | in this case. 

An important difference in the avalanche structures be- 
tween the uniform and disordered Cs systems is the lack 
of periodicity (irregular) in emitter lead voltage and that 
the basins no longer evolve by quenching other basins 
(they overlap). 
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characterizing the scahng of mean linear size and mean 
Rg with area for DC arrays also agree to within errors 
with exponent values from UC arrays. Thus avalanches 
remain essentially compact with a dimensionality of |. 

For LxL^ DC arrays there is no change in the values 
of the exponent a, though the constraining equations 
change to a + 6 = 2.47 and h- ac = 

In this subsection, we have used finite-size analysis and 
been able to successfully relate several finite-size expo- 
nents via scaling relations. Table ^provides a quick sum- 
mary of the values and the context in which they are 
used. Taken along with the fact that these exponents and 
scaling relations help characterize the transition from an 
insulating to a conducting state (the conducting state is 
yet to be discussed), it is reasonable to view the transition 
as a critical transition with associated critical exponents 
and behavior. 
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FIG. 10: Scaling collapse for the distribution of avalanche 
size for UC arrays of size L x L'^^^ . The exponent a and c are 
the same fo r symm etrical systems to within errors, while b 
differs. Fig. |10(b)| plots the data collapse for the mean linear 
size of avalanches with sizes between A and A-l-dA for systems 
of size L xL^''^ with uniform Cs- Exponents are the same as 
those for symmetrical systems. 



Avalanches in DC arrays are not periodic in emitter 
lead voltage and basins don't typically evolve by quench- 
ing other basins - they tend to overlap. This behavior 
is different to UC arrays. By using the three methods 
discussed earlier we find that capacitance disorder does 

shows 



not affect the structure of the avalanches. Fig. 9(a) 
the scaling collapse for the distribution of avalanche sizes 
N(A), for square arrays with nonuniform Cs- The con- 
straining equations in this case are now a+b — 2.8 (as the 
sum of the product of N(A)A for all avalanches ~ L'^^^) 
and b — ac = 2.0. In spite of the presence of capaci- 
tance disorder, the value for the exponent a is the same 
to within errors for the value for uniform Cs . Exponents 



B. Interface motion 

The maximum advance of charge into the system at a 
given Vl can be used to define an interface. Properties 
of the interface can be used to understand other proper- 
ties of the system like fiuctuations in Vr- Some details 
are required about the way we define the interface. At 
a given Vl, there will be some dots onto which electrons 
have not tunneled yet, defined relative to the original 
stable configuration reached by relaxing an original con- 
figuration with < Vi < e/Cs. We refer to such dots 
as zero excess dots. We can define the interface as either 
of the following: (i) contour of leftmost sites along each 
row which has not had an electron tunnel onto it, or (ii) 
the contour of last sites along each row which has had 
an electron tunnel onto it. The two although seemingly 
similar are different in the sense that the second defi- 
nition considers the case where there may be "bubbles" 
of zero electron dots enclosed behind the interface. The 
difference, however, is not significant as the long wave- 
length properties of the interface (e.g., roughness) do not 
seem to depend upon which definition is used. As Vl 
is increased, electrons tunnel onto arrays, via avalanches 
and if electrons tunnel onto a zero excess dots, the inter- 
face advances. The motion of the interface in response 
to a driving force, can be described in the language of an 
elastic medium driven through a random potential. We 
will argue that some quantitative correspondences exist 
in fact. The dynamics of such elastic interfaces through 
quenched disorder has been extensively studied in recent 
years"^^, e.g., CDW, fiux lines in type II superconductors 
etc, fiuid fiow through a porous medium to name some, 
fiux front in thin films of type II superconductors'^^, com- 
bustion of paperiSi^. 

We define the roughness (width) of the interface as the 
square root of the mean of the square of the fiuctuation 
from the mean position. On increasing Vl the interface 
advances further into the system and gets rougher. As 
charge builds up behind the interface, the advance of the 



12 



0.4 
0.35 

0.3 
0.25 

0.2 
0.15 

0.1 
0.05 




* W=16 

- W=25 - 

- W=36 
□ w=49 
■ W=64 



0.05 0.1 



0.15 0.2 
tAV^ 

(a) 



0.25 0.3 0.35 



0.4 
0.35 

0.3 
0.25 

0.2 
0.15 

0.1 
0.05 




* W=16 

- W=25 - 

• W=36 

□ w=49 - 

■ W=64 



0.2 



0.25 



0.05 0.1 0.15 
tAV^ 

(b) 



FIG. 11: Collapse of the roughness of the interface for L x 
L^^^ systems using the scaling form given by Eqn. 1151 Here 
W is the system width and t is measured by the distance 
of t he mea n location of the interface from the emitter lead. 
Fig. |ll(a)| is for systems with uniform Ce , where the values 
of t he expo nents used in the collapse are a= 0.5 and z=1.5. 
Fig. |ll(by| is for systems with non- uniform Ce and the values 
of the exponents used in the collapse are q= 0.5 and z=1.5 
as well. 



interface is analogous to the growth of a surface due to 
deposition of a material. It is weU known, that such 
surfaces become increasingly rough with time, gradually 
reaching a saturation width. For QDA as the advance of 
the interface is governed by Vl; it plays the role of time, 
which upto a constant factor is the same as the mean 
position of the interface. Using the well known scaling 
forn»2i: 

w{L,t) = L^'Wit/L'), (15) 
we were able to collapse data on the width of the in- 



terface with time onto a single scaling curve Fig. 11(a) 



properties of the interface. Due to the large values of the 
dynamic exponent z (1.5), we were able to study only 
small system sizes with interfaces with saturated width. 
Consequently in order to study steady state properties of 
larger interface lengths, it is prudent to study non-square 
systems like LxLs, thereby permitting a more accurate 
determination of the exponents and hence the universal- 
ity class the interface growth pr ocess b elongs to. 

From the the collapse in Fig. ll(a)| we find values of 
the roughness exponent a = 0.5 and dynamic exponent 
z = 1.5 - therefore the growth exponent (3g — 0.33. This 
is consistent with the roughening of the interface being 
in the KPZ universality class^, where a= ^and z— |. 
The KPZ universality class is consistent with the sym- 
metries of the system, viz., rotation is a symmetry on 
large scalesiS, interactions are short range and the speed 
of the interface advance lacks large fluctuations. In light 
of the assumed lack of spatial correlation of the underly- 
ing charge disorder for dot arrays (statistically Galilean 
invariant) and the fact that the interface will move for- 
ward when the emitter lead voltage is increased by one 
unit (^), a description of the interface advance in terms 
of thermal KPZ equation seems valid. 

Some avalanches involve electrons hoping onto a zero 
excess dot - a new-site. When an avalanche involves 
new-sites, the interface is reconfigured; the distribution 
of the avalanches that involve new-sites provides infor- 
mation on the reconfiguration (advance) of the interface. 
The scaling collapse for the distribution of avalanches 
that have betwee n s an d s-|-ds new-sites for UC arrays 
is shown in Fig. (12(a) I. We define N{s) analogous to 
N{A), where s is the number of new-sites visited in an 
avalanche. For LxLs samples the sum of the product 
N{s)s for all avalanches, scales as the number of dots in 
the array (Ls ). Thus the constraint on the exponents /i 
and S in the scaling ansatz: 



Nis) - L'fjis/Lf^), 



(16) 



We initially used symmetrical L x L systems to study the 



is given hy fi + 6 ~ 1.67. Another constraint is de- 
termined by the scaling of the number distribution of 
avalanches with the number of new-sites for a given bin 
([s, s-|-ds]), with system length as La , which results 
in only one independent exponent in the scaling ansatz. 
Hence S ~ iiv = 0.67. We find that the exponent values 
from the collapse consistent with these constraints. We 
interpret the value of the exponent fi — 0.67 as giving 
the typical number of new sites involved in an avalanche 
of linear length I as 1"-^^. We know that the width of 

2 

the an avalanche of linear length I, is also I3 , which indi- 
cates that the avalanche typically involves one new dot 
for each dot along the width. Thus the interface ad- 
vances smoothly on the average by 1 dot along the width 
of the basin of activity. Fig. |12(b) shows the configura- 
tion of the interface at a given Vl and an avalanche that 
crosses the interface with the portion to the right of the 
interface being the new-sites involved in the avalanche. 
These new-sites will determine the new configuration of 
the interface after the avalanche. 



13 



1000 
100 
10 

1 

0.1 
0.01 



0.001 

















H= 0.67 


512x64 


5= 0.97 




343x49 


v= 0.3 




216x36 • 






- 125x25 = 






64x16 -■ 






27x9 o 







0.01 



0.1 



10 



1000 



100 



10 



0.1 



0.01 
0.001 



















11= 


1.0 


1 


512x64 




0.67 




343x49 « 


t)= 


0.05 




216x36 • 








" 125x25 = 








64x16 ■ 








27x9 ° 









0.01 



(s/L'') 



0.1 

(s/L^') 



10 



(a) 



(a) 





(b) 



(b) 



FIG. 12: Scaling collapse for the distribution of the num- 
ber of avalanches involving between s and s+ds new-sites for 
uniform Ce- The exponents are related by /x -I- 5 = 1.67 and 
5 — /ii) = 0.67. The value of ^ = 0.67 indicates that when 
the interfa ce moves it typically advances by 1 dot spacing. 
Fig. |12(a)| shows an avalanche crossing the interface for uni- 
form c'e"^ Notice that the amount by which the avalanche 
overshoots the interface is of the order one. 



FIG. 13: Scaling collapse for the distribution of the num- 
ber of avalanches involving between s and s-|-ds new-sites for 
L xL^^^ systems with non- uniform Cs- The exponents are 
related hy fi + 5 — 1.67 and 5 — nv — 0.67. The value of 
/i = 1.0 indicates that when a seg ment of the interface moves 
it typically advances by l^^^ . Fig. |13(b)| shows an avalanche 
crossing the interface for non- uniform Ce systems. Typically 
the overshoot of the l'^^^ portion of the interface is l^^^. 



Further information on the movement of the interface 
can be obtained by studying the voltages (V^) at which 
an avalanche that involves new-sites occurs, or equiva- 
lently when the interface advances. We can define AVx 
as the difference in Vl between two avalanches that man- 
age to cross the interface (there may be several avalanches 
that do not cross the interface between two interface 
crossing avalanches). Based upon the assumption that 
the advance within basins should be independent, it can 
be shownffi that AVx scales as W/I3 . 

We've seen how the structure of the avalanches is the 
same irrespective of the presence or absence of disorder in 
the capacitance, even though there are changes of major 
significance i n the m otion of the interface. If however, as 
shown in Fig. 13(a) we attempt a scaling collapse for the 



number of new sites covered in an avalanche the exponent 
values are different from the uniform Cs exponent values. 
The exponent /i has a value 1.0 to within errors. This 
value can be interpreted as follows: the dimensionality of 
avalanches is |, which means for linear size I the width 
is typically Zi. When an interface crossing avalanche 
occurs, the average amount by which it overshoots the 
interface scales a.s , hence covering l^xls new sites. 

This can be seen in Fig. 13(b)| which represents a typ- 
ical interface crossing avalanche in a sample with dis- 
ordered Cs, where the avalanche overshoots the inter- 
face by a significant amount compared to the uniform 
Cs (where the overshoot was of order one spacing) . For 
DC arrays avalanches do not occur with any fixed reg- 
ularity - either spatial or temporal - so a large number 
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of avalanches may occur which do not reconfigure the 
interface, followed by an avalanche that rearranges the 
interface by a large amount. Compared to the smooth 
motion of the interface in arrays with uniform Cs, the 
motion of the interface for DC arrays is rather jerky. It 
is important to mention that the motion appears jerky 
locally, but at a coarse-grained scale and on average the 
velocity of the interface is well defined and smooth till it 
reaches the collector lead. 



C. Threshold voltages revisited 

Similar to one-dimensional systems the Vt for two- 
dimensional systems scales linearly with the system 
length. For two-dimensions however, there is an addi- 
tional dependence on the width of the samples, which 
can be understood using the concepts of basins and in- 
terface advance from earlier subsections. It also helps 
explain voltage fluctuations. 

With increasing Vl the interface advances further 
along into the system until fi nally electrons reach the col- 
lector lead at Vt- Fig. 14(a) shows how Vr normalized by 
system length (L) depends upon the width of the system 
studied. When Ni,{L) > 1, in addition to fluctuations 
within a single basin, Vt is determined by the basin that 
moves the interface to the right lead the earliest. With 
increasing Nb{L), the expectation value of the maximum 
advance of the interface at a given Vl increases, explain- 



ing the decreasing value of The sample-to-sample 
fluctuations in Vt can be attributed to the roughness of 
the advancing interface. We saw that the roughness ex- 
ponent a for the interface was ^. Assuming a value of 

z= |, for L X W samples, where W = the interface 
reaches its saturation width given by W^, which is . 
As shown in Table UTI fluctuations in Vt agree with this 
picture. Local values of t he thr esho ld fluc tuation expo- 
nents are plotted in Fig. 15(a) and |15(b) and they are 
consistent with the scenario depicted. 

A finite-size scaling length can be defined in terms of 
the characteristic fluctuations in Vt- a{VT)/VT can be 
thought of as defining the scale characterized by L"^/"^^, 
where vt is the exponent characterizing the finite-size 
effects on the transition and for both uniform and disor- 
dered Cs systems has a value of |. 

Analogous to the ID systems, we investigated the re- 
sponse to changed boundary conditions - measured by 
the change in right lead voltage Vr -hy measuring differ- 
ence in Vt, as a method of probing the disorder energy 
scale. It is shown''" that {AVt)avtt^o scales as La, i.e., 
for 2D when Vr changes, it changes on average by a value 
given by La. We have discussed the mapping between 
Eden growth (which is in the KPZ universality class) 
and DPRM2&. Using the analogy, the maximum point of 
advance of the interface in our systems can be mapped 
to the ground state of a DPRM. It is known that the free 
energy fluctuations of the ground states, both for sample- 
to-sample fluctuations and higher order excitations scale 
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FIG. 14: The dependence of Vt on the numb er of b asins, 
N5(L) for uniform Ce systems is plotted in Fig. |14(a)| The 
dependence of Vr on the number of basins (Ni,(L)) for non- 
uniform Ce systems is plotted in Fig. |14(b]| 



as La. For systems whose ground state (maximum ad- 
vance of interface) is unable to overcome the increased 
voltage of the right lead (a change in boundary condi- 
tions) the next lowest energy state (interface position) 
needs to be enough to overcome the changed boundary 
condition; the energy of which is typically La higher than 
the ground state. The La behavior can be understood 
without invoking the mapping between Eden growth and 
DPRM. We saw in the subsection on interfaces, that the 
mean voltage increment to move the interface so as to 
have a new maximum position scaled as La . When 
the right boundary condition is changed, either the last 
avalanche is able to overcome the increased right lead 
voltage (in which case AVl = 0) or requires an increase 
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FIG. 15: EflFective exponents for the fluctuation of threshold 
voltage for arrays with uniform capacit ance an d array s with 
disordered capacitance are plotted in Fig. |15(a)| and Fig. |15(b)1 
respectively. Sample-to-sample fluctuations in Vr for nonuni- 
form Ce arrays are similar to the uniform Ce and scale as 
L" "^'^ to within errors for all values of Nb{L). 



in Vl, in order to surmount the barrier at the right lead. 

In addition to a well defined critical point, a true con- 
tinuous phase transition is characterized by the fact that 
there aren't any characteristic length scales in the sys- 
tem, i.e., fluctuations take place at all length scales. Ob- 
viously this is not true for finite size systems - where 
possibly all length scales upto the system size, but no 
larger can be present. This system-size dependent cut- 
off explains why we see finite-size deviation from the true 
(infinite-size) values. There is a systematic dependence 
of these deviations with the system sizes studied. By 
examining these systematic dependences on the scaling 
exponents, we try to extrapolate to the infinite-size value 
of the exponents. The fact that a system-size indepen- 
dent behavior is possible over a range of system sizes 
(e.g., the collapse of several system sizes onto a single 



TABLE I; Symbols used and comparison of values for uniform 
and non-uniform Ce samples. 

Symbol Used in uniform Ce non-uniform Ce 

^~c N(A) vs A 1.7, 1.3, -0.43 1.7, 1.1, -0.55 

K 1 vs A 1.0, 0.63, 0.58 1.0, 0.67, 0.67 

(3g, z Family- Vicsek scaling 0.5, 0.33, 1.5 0.5, 0.33, 1.5 

11, 5, V N(s) vs s 0.67, 1.0, -0.3 1.0, 0.67, -0.05 

9 fluctuations in Vt 0.33 0.33 

T mean of nonzero AVt 0.0 0.3 



a, b 
p, a 
a 



TABLE II: Fit values for for uniform and non-uniform Ce 
samples, with different number of basins (Nt). For (t(Vt) = 
AL^ both A and 9 were fit parameters. 

Nb{L) uniform Ce non- uniform Ce 



0.33 ± 0.01 
0.34 ± 0.01 
0.35 ± 0.01 
0.36 ± 0.02 



0.33 
0.34 
0.34 
0.35 



0.01 
0.01 
0.01 
0.01 



curve) is the crucial indication that finite-size scaling ap- 
proach is successful, which in turn is an indication of a 
phase transition. Hence the claim that Vr can be viewed 
as the critical point of a phase transition. 



IV. 2D ARRAYS AT THRESHOLD 

Vt is defined as the lowest lead voltage at which there 
exists at least one dot in the column adjacent to the 
emitter lead, onto which electrons can tunnel and subse- 
quently find a way onto the collector lead. Having estab- 
lished the existence of a threshold voltage in the previ- 
ous section, our aim in this section is to understand the 
current conduction in the array at exactly the threshold 
voltage. A few samples of the first conducting path at 
Vt are shown in Fig. 1161 from which it can be seen that 
there are frequent splittings and possible recombination 
of paths, leading to an overall complicated geometry and 
topology. We investigate the structure and the current 
carrying capacity of the ground state paths. We will find 
immediate application of our understanding in the next 
section, when we investigate the I-V characteristics of 2D 
arrays. In the next subsection we present the details of 
how the first conducting path at Vt is determined. We 
then present results on the transverse deviations (mean- 
dering) of the path - characterized by a wandering expo- 
nent (C) - followed by a discussion of the main structural 
features of paths. Finally, we discuss the current density 
profiles at the end-points and establish a connection be- 
tween structure and current density profiles and compare 
the ground state path for UC arrays and DC arrays. 



A. Computing the ground state path 

To describe the ground state path, in addition to de- 
termining the location of where electrons fiow, we are 
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(c)c 

FIG. 16: A few ground state path samples illustrating the 
merges and splits along a single path. Some splits are "effec- 
tive" in that paths do not recombine, for example, the middle 
split in figure a. Some paths have multiple splits but none 
go onto persist till the end (figure c). Then there are paths 
that split early on and then wrap around - due to periodic 
boundary conditions, and merge with the original path (fig- 
ure b). This results in a mouth width equal to the system 
width. The gray-scale is an encoding for the current density 
at a given dot. 

interested in determining the relative probability of an 
electron tunneling through a given site. We will use rel- 
ative probabilities as measured by current densities (j) 
[to be defined in Eqn. (|17|l ] and not macroscopic currents 
(I). We use a transfer- matrix approach to determine the 
relative probability of electrons tunneling through dots. 
This involves computing the probability of being in state 
i, using known probabilities of being in all possible previ- 
ous states j and the transition probability of going from 
the states j to state i. Due to the numerical uniqueness 
of the random potential at each site, there is in practice 
only a single dot onto which electrons can tunnel from 
the emitter lead at Vt- We assign this dot, which is at the 
same potential as the emitter lead, a current density of 
1.0, as all current flowing onto the array passes through 
this head dot. It can be shown that an electron cannot 
tunnel onto any other dot in the left most column other 
than the root of the spanning avalanche'^° - the head dot. 
Thus all other dots in the leftmost column are assigned 
a probability of 0.0 (the boundary condition). 

As electrons can tunnel onto a dot only from dots that 
have a higher electro-static potential. Hence in order 
to compute the current density of a dot (probability of 



an electron flowing onto a dot), the current density of all 
neighboring dots which have a higher potential should be 
known. The current density of any dot i {ji ), is computed 
as the product of the current densities of dots in the 
immediate vicinity of i with the probability of tunneling 
from the neighboring dot onto dot i, summed over all 
dots: 

where jj is the current density of dot j and Pj->i is the 
probability of tunneling from a dot j to i. Pj->i is com- 
puted as, 

p.-. = r,-./rLt, (18) 

where Tj-^i/T''^^^ is the ratio of the tunneling rate from 
dot j onto dot i over the sum of all outgoing rates from 
dot j. As the probability of tunneling onto a dot from 
a dot at lower potential is zero. Thus starting with the 
head dot with a current density of 1.0 and sorting all dots 
in decreasing order of potential, the current density is 
computed for the dot with next highest potential. As the 
current densities and the probabilities of tunneling are 
known for all dots from which electrons can tunnel onto 
it, the current density for the new dot can be determined 
using Ean. lfTTji . 

A special case of ji is jL{i), which is defined as the 
current density from the ith dot in the last column onto 
the collector lead. It is useful to note that the sum of 
the j's along any column can be greater than 1.0 (e.g. 
when there is intra-column hopping), but the sum of all 
current densities between adjacent two columns must be 
equal to 1.0 (essentially a current continuity equation). 
Thus the sum of all jL(i) will be 1.0. 

There isn't a simple connection between the current 
densities computed using our approach and the actual 
macroscopic currents that can be carried by a path. The 
current densities approach taken here, provides informa- 
tion on the relative proportion of the current that would 
tunnel off the dots onto the collector lead, i.e., be carried 
along different paths, but says nothing about the exact 
values corresponding to a given path. It is possible for 
example, that at threshold, a simple non-splitting path 
conduct more current than a highly complex path with 
many splittings and recombinations. 

B. Ground state path properties 

1. Path meandering, widths and energy fluctuations 

The number of end-points (riep) is deflned as the num- 
ber of dots in the last column - adjacent to the collector 
lead - which have a nonzero value of ji (strictly speak- 
ing, a non-zero value of jhii))- As shown in Fig. 1171 it 
becomes exponentially less probable to find paths with 
an an increasing number of end-points. It is relatively 
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FIG. 17: The probability distribution of the riep taken deter- 
mined using approximately 75000 disorder realizations. There 
is a good fit to an exponential line, implying that there is an 
exponentially decreasing probability of a path having higher 
number of end-points. Other system sizes have a similar dis- 
tribution. 



simple to implement a tracking algorithm that by work- 
ing downwards from the head node computes the trajec- 
tory of the electron hopping and determines the number 
of transverse (inter-row) deviations en-route to the end 
nodes. The deviation of end-points is of interest and 
for the i*'* dot is given by The current density 



weighted transverse deviation can be defined as: 

0i = 3LyL{i) 



(19) 



and the current density weighted mean transverse devia- 
tion as: 



(20) 



As there are often more than a single end-point, thus 
the weighted mean $ determines the location of the ef- 
fective end-point of the path and thereby the deviation 
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FIG. 18: Plot of the value of the effective value of the fluctu- 
ation of $ with L which gives the wandering exponent. The 
wandering exponent approaches a value consistent with |. 



of the path from the head node. The value of $ av- 
eraged over many samples will be zero, as there is an 
equal probability of the effective end-point being on ei- 
ther side. A look at the values of <I> over many disorder 
realizations, reveals essentially a Gaussian distribution 
about the mean value and the standard deviation of the 
distribution grows as , where C, is the wandering expo- 
nent and is found to have a value of | . This sets the scale 
for typical sample-to-sample fluctuations of the effective 
end-point. Fig.^|plots the local value of the exponents, 
which in general is a useful way to get a handle on the 
finite-size dependencies of the exponentsi^ 

Computing the transverse deviation associated with 
dots helps calculate the width of the "mouth" of the 
path in the last column, which is defined as the differ- 
ence of the transverse coordinates of the extreme end- 
points. It is of interest to understand how the width of 
the mouth grows with system size. The wandering expo- 
nent provides insight into the typical fluctuations of the 
location of current density weighted end-point but does 
the width of the mouth grow with the same exponent? 
The mean value of mouth- widths for UC arrays is shown 
in Fig. ll9f a). The growth in the width is consistent with 
La - definitely different from La. This is indicative of 
a situation where the location of the effective end-point 
fluctuates as L^ but the distribution of the extreme end- 
points around the effective end-point gets "smaller" rela- 
tive to the effective end-point fluctuations. The increase 
in the mean width of the mouth tells us that paths that 
require the same potential difference as the ground-state 
path to within 0(1), are to be found upto La around 
the effective end-point. The fluctuations in the width of 
the mouth is plotted in Fig. Iiyr bl. The fact that the 
fluctuation in the width scales as La, indicates that the 
extremities of the mouth are determined by randomly 
juggling end points. We've discussed the fluctuations in 
the threshold voltage in the s ection UlI CI we remind the 
reader, that as shown in Fig. 15(a)| and Fig. 



15(b) the 



fluctuations in the threshold voltage scales as La . 

As a quick consistency check that the wandering ex- 
ponent is not dependent on the width, we compared 
the wandering exponent for systems with TV;, = 4 
(LxA^fjL^/^) with those of systems with Nb = 1. Al- 
though there are significant differences at smaller system 
sizes, for larger system sizes the boundary effects become 
less significant. For systems with Nb =4, the convergence 
to the L^^^ is sooner than for Nb =1, indicative that 
boundary effects dominate at small system lengths. 



C. Path geometry and topology: Gap sizes, merge 
lengths and typical splitting distances 

So far our understanding of the structure of the ground 
state path is that there are possibly many branchpoints 
leading to multiple end-points. The current density 
weighted transverse deviation leads to an effective-end- 
point, with sample-to-sample fluctuations of La and 
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these questions, will help understand several important 
length scales characterizing paths at threshold. It will 
also provide additional metrics to compare the ground- 
state paths of systems with different disorders (UC, DC 
and RT). 

In order to compute the size of gaps separating the 
end-points and to compute how frequently and over what 
length scales paths split we need to look beyond the 
transverse deviations of each point. We need to track the 
complete trajectory electrons may take before tunneling 
onto the end-points. This involves knowing the location 
of all branchpoints along the path , where a branchpoint 
maybe defined as a location at which a split occurs, i.e., 
where there is more than one neighboring dot onto which 
electrons can tunnel. There are many splits along the 
path; the majority of the splits along the path do not 
survive and merge a short distance after splitting. In the 
thermodynamic limit, not all branchpoints are of inter- 
est, but only those branchpoints that go on to produce 
end-points - do not merge after splitting. We find the 
last possible location that is common to the trajectories 
associated with the two end-points of interest. This point 
is referred to as the effective branchpoint corresponding 
to the two end-points. Having thus determined the lo- 
cation of the effective branchpoints for all pairs we can 
compute the position at which paths to the end-points 
last overlap. Equivalently, this location can be used to 
establish a lateral distance from the collector lead that 
paths from end-points a transverse distance Ay^ apart 
will most likely merge by. 

We then compute the mean value of the lateral length 



FIG. 19: The width of the mouth of the path is defined as the 
[max{yL{i)) — min{yL{i))]- In Fig. ll9f a'l the left y-axis gives a 
plot of the width as a function of system size; using the right 
y-axis gives the local slope of the width with system size. 
As can be seen the val ue of the exponent settles to around 
0.37 ±0.03. Fig |19(b)| shows the results for the fluctuation 
of the width of the mouth of the path. The growth in the 
fluctuations of the width are consistent with Ls . The left- 
hand y-axis represents the sample-to-sample fluctuations of 
the width, whereas the right-hand y-axis gives the value of 
the local slope. 



mouth- width which scales as L^. As a consequence of 
the many interspersed end-points between the extremal 
end-points of the mouth, the mouth has a fine structure 
not accessible by studying only the transverse fluctua- 
tions of the effective end-point and widths. We would 
like to understand the details of fine structure of the 
mouth, viz., to understand if any two physically contigu- 
ous end-points are part of the same path segment or if 
they belong to two different path segments. In general, 
if two end-points belong to different segments, typically 
how far back did those segments separate? Answers to 
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FIG. 20: Fig. |20(a)| is a plot showing the mean lateral length 
for all gaps of size A. The mean lateral length of a gap of size 
A, can be thought of as the typical distance at which the path 
to the head-node from two end-points a distance A apart will 
have merged. The fit indicates that the lateral length scales 
as 3/2 of the gap size. The plot also shows a comparison of 
the mean lateral lengths for the maximum and minimum gaps 
for an effective branchpoint. 
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of gaps of size A. Given that the wandering exponent has 
a value of |, one would expect that the paths to two end- 
points separated by Ay^ would be typically joined upto 
a distance Ayi,^^'' from the collector leads. This is anal- 
ogous to the typical separation of the optimal paths to 
two ends of a DP RM th at arc Ayj^ apart, viz., Aj/^^^''). 
As shown in Fig. 20(a) our findings are in good agree- 
ment with this expectation; the mean lateral length for 
all the gaps for the n end-points is used in the fit. 

The fact that the path structure of QDA is similar to the 
scale-invariant tree structure of DPRM, is further indi- 
cation that ground state paths of QDA are in the same 
universality class of DPRM. The significance of this con- 
clusion will be explored later. Also plotted are the mean 
lateral lengths of the maximal and minimal sized gaps of 
an effective branchpoint, which scale like | and | respec- 
tively with the transverse size of the gaps. 

The gap between two end-points is defined as the num- 
ber of the intermediate dots separating them. Thus 
for two physically adjacent end-points (irrespective of 
whether they belong to the same path segment or not), 
the gap is defined to be of zero size. We computed the ef- 
fective branchpoints for all pairs of end-points (there are 
"^"^""'"^ pairs for n end-points) and computed the lateral 
and transverse sizes of the gaps. The results of the prob- 
ability distribution for a fixed system size arc shown in 
Fig. 21(a) and Fig. 21(b) The probability distributions 



represent the simple fact that large gaps resulting from 
earlier permanent splittings of the path at the threshold 
become less probable. As the path at threshold and the 
end-points have become reachable within the last 0(1) 
increase in potential, all path segments at threshold must 
be equal to each other to within voltage of 0(1). Con- 
sequently they will overlap for the most part. We have 
seen that the sample to sample fluctuation in the thresh- 
old voltages scales as La which should set the scale for the 
typical difference between non-overlapping paths, thus at 
threshold we expect typically ^^t^ fraction of paths seg- 

ments to not overlap. Consequently if for a given system 
size, we were to plot the mean of all lateral sizes of the 
gaps as a fraction of the system length, we'd expect to 
see a L~3 dependence. 

By studying the distribution of the location of the ef- 
fective branchpoints, and found that it became increas- 
ingly improbable that they would be located closer to the 
emitter lead^. It is also useful to compute the number of 
effective branchpoints (depth) encountered on a path to 
an end-point and how the depth varies for the different 
end-points in a given sample. This tells us if the paths to 
the end-points are typically similar, as well as permitting 
us to determine a correlation between physical proximity 
of end-points and difference in depths. Thus we deter- 
mine the average number of end-points for a given mean 
depth of end-points. If the tree structure was perfectly 
random the number of end-points would grow like the 
square root of the depth. On the other hand for an essen- 
tially unbalanced trees (where all the splittings take place 
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FIG. 21: Fits to the probability distributions of the gap sizes 
(Fig.EHa)) and the lateral length of the gaps (Fig.EHb)) for 
a single system size. Although not shown, there is not a sys- 
tem size dependence at which probability decreases (though 
there is a system size dependent cut-off). Fits indicate the 
probability of occurrence of a gap of size A decreases as A~^. 
A similar dependence is observed for I, which is consistent 
with expected probability distribution for I computed using a 
variable transformation from A to / (where A ~ l^^^). Given 
the distribution of I, the chance that paths that spht, will sur- 
vive as independent paths all the way to the end gets smaller 
the earlier they split. 



on the path to one particular end-point), the number of 
end-points, n^p grows linearly with the mean depth. For 
an essentially balanced tree each path splits essentially 
with equal probability in which case the number of end- 
points grows as some number to the power of the mean 
depth (for a perfectly balanced tree this would be the 
number would be 2). As shown in Fig. [221 we find that 
the mean depth grows logarithmically with the number of 
end-points, characteristic of an essentially balanced tree. 

Finally we'd like to determine if the trees (representing 
the ground-state paths) are spatially homogenous and if 
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to the nodes typically have a net difference of one addi- 
tional branchpoint with every scale of two increase in the 
separation between nodes. 

It is fair to assume that path segments that are not im- 
mediate descendents of the same parent are essentially in- 
dependent; every pair of end-points with a Adcg > 1 can 
be considered independent and thus the plot in Fig. 1231 
provides a measure of the number of independent paths 
segments (channels) that reach the collector lead. This 
could possibly be experimentally verified by studying the 
spectrum of the discrete current at the collector lead. Us- 
ing this definition of 'independent' paths, we find that 
the number of independent channels increases logarith- 
mically in the transverse direction upto the width of the 
mouth. 



FIG. 22: In Fig. 1221 the mean value of the dcg of the end- 
points is plotted as a function of riep. As the number of end- 
points increases the mean depth of the end-points increases 
logarithmically. This is indicative of an essentially balanced 
tree. For a perfectly balanced tree the number of end-points 
increases exponentially with the depth, whilst for an hierar- 
chically split tree (in the limit of a perfectly unbalanced tree) 
the mean value of the depth increases linearly with the num- 
ber of end-points. 
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FIG. 23; Aj/L is the transverse separation of two end-points 
on the path and Adcg is the difference in the dcg between the 
nodes. Given the logarithmic increase in the value of Adcg 
with transverse separation, the paths to the nodes typically 
have a net difference of one additional branchpoint with every 
scale of two increase in the separation between nodes. 



the path segments from the different end-points to the 
head node, are essentially similar in the number of effec- 
tive branchpoints encountered. To do so, we investigate 
the difference in the depth (on the path to the end-point), 
represented as Adcg, of end-points separated by a trans- 
verse distance Ay^. As shown in Fig. 123 the difference 
in depths increases logarithmically with transverse sepa- 
ration of end-points. Given the logarithmic dependence 
of Adcg on transverse separation of end-points, the paths 



D. Current densities 

We have found that the structure and the topology 
of the first conducting path to have several interesting 
features and characteristic length scales. An important 
question is what is the profile of the current densities 
at the end-points? Also, what does the fluctuations of 
the current densities within the mouth tell us about the 
overall structure of the paths? 

To address what the current density values for a pair 
of end-points tell us, we take and as the values of 
the larger and smaller current density (for the two end- 
points in consideration) respectively, and as a measure of 
the difference in the number of splittings encountered for 
two end-points define An^ as: 

An, = log(^) (21) 

In Fig. we plot the value of An, as the transverse 
separation between end-points increases. The best fit to 
the largest system size considered (L=3375, W=225) is 
consistent with a logarithmic dependence over the range 
1 to about 20. From the logarithmic increase in Aris with 
transverse distance (Ay^), it follows that j> = Ay^j^, 
i.e., with increasing distance between the two points con- 
sidered, the larger current density (j^) tends to get larger 
relative to the smaller (j^) - increasing linearly in Aj/^. 
We know from insight gained from the structure of the 
paths, that as the transverse distance separating two end- 
points increases, the paths taken to the two end-points 
separate earlier. If the paths to the end-points after sep- 
aration typically undergo the same number of current 
splittings, then on average there would not be any varia- 
tion in Arig with distance; but given the slow but definite 
distance dependence, it is consistent to conclude that one 
path undergoes more current splits than the other, and 
that for end-points separated by a greater transverse dis- 
tance, the correlation in current densities will be less than 
for those end-points which have greater overlap in their 
paths. 
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FIG. 24: Aris is a measure of the difference in the number 
of s plits in the paths traversed to get to two end-points. In 
Fig. |24(a)| the dependence of the mean value of An^ on the 
separation Ayi, is plotted. As can be seen initially there is a 
logarithmic dependence of An^ on Aj/_t (the fit shown here is 
for the largest system size) before gradually crossing over to 
a Aj/L independent value. Data for widely differing system 
sizes (from L=343 to L=3375). To within errors, the plateau 
value appears to be independent of the system size. 



It is useful to point out the similarity betweenthe log- 
arithmic dependence of An^ on Ay^ as in Fig. 24(a) and 



the logarithmic dependence between the dcg on Ay^ as 
shown Fig. 1231 In general, given the similarity in the 
properties of the current densities at the end-points and 
the structure of the path, it appears to be the case that 
the effective branchpoints not only determine the struc- 
ture of the paths but also play a role in determining the 
current profiles of the end-points. 
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FIG. 25: The values of the local slopes for the wandering ex- 
ponent compu ted for arrays wit h disor dered capacitance are 
plotted in Fig. |25(a)| From Fig. |25(b)| deg increases logarith- 
mically with Hep. Both are essentially indistinguishable from 
UC arrays. 



E. QDAs with capacitance disorder 

In this subsection, we study the properties of the first 
path for DC arrays and begin by investigating the trans- 
verse deviations of the path and the structure of the 
mouth for ground-state paths. As shown in Fig. |23fa), 
the wandering exponent gradually approaches the value 
of ^ = I for larger systems, which is similar to UC ar- 



rays. In Fig. 25(b) the relationship between dcg and 



is shown to be logarithmic (recall that the probability of 
occurrence decreases exponentially as Uep increases). 

In Figs.|2niwe plot the distribution of the gaps and the 
mean lateral distance of splitting for gaps of size A. Both 
the distribution of the gaps sizes (and thereby lateral size 
of the gaps) and the mean lateral distance dependence on 
gap sizes are similar to UC arrays. 

In addition to the structure of the path, current flow 



properties are also indistinguishable to UC systems as 
shown by the sample averaged fluctuations of the current- 
density weighted transverse location^. From the data 
as presented in this section, ground state paths are effec- 
tively indistinguishable from the UC. It is highly unlikely 
that any further investigations will indicate any signifi- 
cant differences between the ground path structures for 
the UC and DC systems. 

DPRM is controlled by a zero fixed point thus the 
ground state (lowest energy) strongly determines the 
properties of the system. Given the fact that the ground- 
state conducting path is in the same universality class as 
the DPRM, one would expect excited conducting paths - 
those with energy higher than the ground state at higher 
voltages as well as ground-states at non-zero tempera- 
tures - to be strongly infiuenced by the structure and 
energetics of the first conducting path. Thus the connec- 
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FIG. 26: Comp aring t he gap distribution and spacing for UC 
and DC in Fig. |26(a)1 and Fig. l26f b'). indicates that they are 
essentially indistinguishable. 



tion between QDA and DPRM in addition to providing 
an idea of the structure of the ground state path at zero 
temperature, indicates that an extension of the approach 
used here to study the ground state path might possibly 
be used in determining sensitivity to boundary conditions 
and temperature changes of the ground state paths. The 
latter is of significant practical importance. Given the 
putative similarity between QDA and DPRM we can use 
results obtained in the DPRM case to predict a temper- 
ature sensitivity: namely that the ground state configu- 
ration is sensitive to temperature changes and will most 
likely rearrange. As to whether this is sufficient to change 
any scaling properties will require explicit numerical and 
analytic work. 



V. CONDUCTION IN 2D ARRAYS 

In the previous sections, we saw how the threshold 
voltage can be viewed as the critical point of a continu- 
ous phase transition and explored the associated critical 
phenomenon at and below the critical point. This sets the 
stage to address the next, and arguably most important 
question in our investigation of disordered QDA - the 
nature of the critical phenomenon for voltages above the 
critical point. Based based on the strength of the driving 
force {V) relative to the strength of internal interactions 
and disorder strength, roughly three distinct regimes can 
be identified. The first regime can be thought of as when 
the scales of disorder, interaction and driving force are 
all similar. In this regime the role of disorder is generally 
crucial and the interactions between the many degrees- 
of-freedom result in strong deviations from a mean-field 
behavior. This regime typically occurs when V is very 
close to the threshold voltage. A second regime lies at 
the other end of the spectrum, where the driving force 
is extremely strong compared to the strength of disorder 
and interactions between the degrees of freedom; in this 
regime the disorder and interactions become irrelevant 
and the system is driven into a linear response mode. Our 
primary focus will be on the investigation of the critical 
behavior and dynamical response close to the transition 
- corresponding to the first regime. We will study the 
dynamic response by computing the I-V characteristics 
for a range of different systems sizes. Details of theory 
and implementation of our numerical simulations can be 
found in Ref. [30l |. 

The relative strength of the interactions and disorder 
in turn has been used to broadly classify two widely dif- 
fering types of collective transport: weak disorder rela- 
tive to the strength of interactions most likely leads to 
an elastic structure without breaking up; an example of 
which are CDW. In general when the disorder is strong 
relative to the interactions, the elastic structure breaks 
and transport is far more inhomogeneous and plastic like. 
Examples of transport in such a plastic regime include, 
the flow of a non- wetting fluid in porous medium.'^'', the 
transport of strongly pinned two-dimensional Abrikosov 
flux arrajiS^, driven collective transport of neutral carri- 
ers in randomly varying trap a^&i?? and the flow of a fluid 
with no elastic interactions flowing down a rough inclined 
plane"*" (the dirty windshield problem) . We will find that 
conduction in the low ly regime is plastic-like, i.e., along 
well defined narrow channels. 

Recall that below Vr, the concept of an advancing 
elastic interface - defined as the contour of maximum 
advance of charge along a given row was useful. As a 
consequence of our definition, this elastic interface is no 
longer well defined at driving voltages above threshold, 
and thus not the interface that tears and results in plas- 
tic flow. This leads to an interesting situation where 
as a consequence of asymmetry around threshold, the 
variables and description of the system on the opposite 
sides of the critical point are different; consequently the 
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same exponents are not valid both above and below the 
transition point. This is unlike many continuous phase 
transitions, especially equilibrium (e.g., two-dimensional 
Ising magnets in the absence of an external field) but even 
non-equilibrium phase transitions (e.g., CDW) where the 
same exponents with possibly the same values character- 
ize the critical regimes on either side of the critical point. 

It is instructive to review the MW scaling hypothesis 
originally presented in Ref. [l^ to understand current 
flow in a two-dimensional arrays before discussing the nu- 
merical results. Similar to one-dimensional arrays, any 
current-carrying channel at a given emitter lead voltage 
Vl, there will be ^^iK^ extra charges on average. The 

locations of these extra charges can be viewed as charge 
steps relative to the threshold configuration. Where ex- 
actly these extra charges are located on the channel de- 
pends on the underlying disorder; typically charge down- 
steps are where the tunneling rates are sufficiently smaller 
than the mean tunneling rates. The location of these 
steps give the most likely locations of a split in the path; 
and thus can be used to define a correlation length ^||, 
where = ■jyj-zr^pjQ^ , where L is the linear dimension 
between the emitter and collector leads. We have seen 
that the transverse deviation {£^±) of a path segment of 

length is given by ^jj^'^- Also Vr L, thus ^|| ~ 

and therefore ~ sets the scale for separation 

between channels before splitting. The number of chan- 
nels {Nch) at the collector lead will thus be given by 
where W is the width of the array. Under the assump- 
tion that each channel reaching the collector acts as an 
independent one-dimensional current-carrying chain, the 
current in a channel is Ich ~ i^- Thus the total current 
carried by the array will be given as: 

/ ~ iV,^ X I,h ~ J^'/' (22) 

It is important to discuss some of the assumptions that 
the MW hypothesis depends critically on. 

Firstly, that each channel behaves like a one- 
dimensional array and the current in the ID array grows 
linearly with ly. Secondly, that number of channels grows 
like ~ v^^^ , which in turn is dependent upon two assump- 
tions. The first is that the transverse deviations grows 
like P^^, where I is a linear dimension of the path. We 
have extensively verified this to be true at Vr', it is fair 
to assume that it is above Vr too. The second assump- 
tion is that the most likely effective splits - splits that do 
not result in a recombination - take place at the charge 
down-steps. This has been harder to verify rigorously; at 
best we find for arrays at Vr that the sample-averaged 
probability of effective splits decreased as j. We find 
that the fiuctuations in the current carrying capacity of 
one-dimensional arrays decreases as ~ 

It was originally predictedi^ that to observe the true 
exponent arrays larger than 400 x 400 would be required. 
We will go onto show that the linear dimension required 
before the "true exponent" might be observed is probably 
an order of magnitude larger than initially estimated. A 



significant portion of the remainder of this section will 
be devoted in support of this statement. 



A. Simulation results, analysis and discussion 

From our analysis of the avalanches and path mean- 
dering, we developed an understanding that on average, 
a LxL3 sized system contains one independent basin and 
thus one channel. Before we can understand the current 
carrying capacity of several interacting channels, it is im- 
portant to determine how the current carrying properties 
of a single channel changes with driving voltage. Conse- 
quently, we will most often investigate the I-V curves of 
asymmetrical systems of length L and width Ls . Due to 
finite size effects and computational cost, this also hap- 
pens to be a practical approach. 

We are investigating the scaling properties of a hy- 
pothesized power law between I and v, thus rather than 
perform a coarse grained fit to the entire I-V and gener- 
ate a single value for the scaling exponent, we compute 
"local" values of the exponent /? as a function of i^. We 
find this a more useful and meaningful representation for 
determining how the scaling relation between the I and 
v changes. The procedure although helpful, is not suffi- 
ciently sophisticated to be a complete replacement for a 
rigorous fitting, as error bounds with confidence intervals 
are not easily determinable from this approach. 

The local exponents for UC arrays are plotted in 
Fig. I27r b). from which there is a clear dependence on 
system size for the local exponents. There isn't a range 
of v, however small, where the local exponent curves for 
all the different system sizes lie on a single curve, as would 
be expected for a valid scaling regime. Thus it is difficult 
to claim that there is a unique single value of the local 
exponent for all sizes, even over the smallest regime of v. 
At lower values of the statistical noise starts to dom- 
inate and the true value of the local exponent becomes 
unclear. 

The aim of rigorously verifying the MW scaling hy- 
pothesis numerically does not appear to be easily at- 
tainable with available computational resources at the 
present moment; thus it remains open, as to whether the 
MW scaling hypothesis is valid. If we assume that MW 
is the correct hypothesis, we can at best determine the 
constraints on system sizes and values of to establish a 
regime for the validity of the hypothesis. This is some- 
what analogous to determining an upper bound of the 
reduced variable upto which critical behavior can be ob- 
served: 

It is known that for CDW one has to be within 
/ < 0.01— (where / — -^^-^r^) of the critical point in 
order to observe associated critical phenomenon. Simi- 
larly for high-Tf, superconductors (copper oxide) in three- 
dimensions, by some estimates^ the critical region exists 
for t — 10^^ where t = C^^'^") The quoted estimates are 
from analytical calculations and supported by numeri- 
cal data. With the caveat that it is much harder to 
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FIG. 27: I-V curves for 2D arrays with offset charge disor- 
der for a range of system sizes. The Piocai for the I-V curves 
in Fig. I27f a') is plotted in Fig. I27f b'). For a true power law 
scaling, the value of Piocai for different system sizes should 
overlap. As can be seen this does not happen for values less 
than V = 0.1. Following the MW scaling hypothesis, we ex- 
pect a plateau at values less than u = 0.1. The inability to 
see clearly a definitive plateau is primarily due to large finite 
size effects. 



estimate correctly using numerical data alone, we haz- 
ard an estimate of the critical region for QDA solely on 
numerical data. From the plot of the local slopes for 
the largest two-dimensional QDA simulated [L = 2744 
and = 196 in Fig. EHb)) and the plot of Piocai for 
largest one-dimensional systems (L= 2000) a similar up- 
per bound would be somewhere between v = 0.01 — 0.1 
for LxLa arrays. As can be seen from Fig. I27r b). there 
are strong indications of a plateau in Piocai , albeit over a 
small region - for values of v around 0.1 - and only for the 
largest arrays. In addition, from the very brief flattening 
out of (3iocai for 1000 X 100 arrays around i/ = 0.1 before 



dipping, it is conceivable that for values of v < O.l, the 
"true exponent" value lie somewhere in between 1.5-2.0; 
this is consistent with the hypothesized value of Piocai 
= 5/3. If at all, this will be the critical region and the 
likely value of the critical exponent. 

It is clear that simulations of even larger system sizes 
will be required to observe a plateau for at least a decade 
'm V - which is the really the minimum range over which 
a power law should be observed before definitive claims 
of scaling are valid. As mentioned simulations of systems 
large than 2744 x 196 are currently computationally not 
feasible. 

For values of > 0.1 the values of (3iocai are influenced 
by a crossover to a peak value of approximately 2.0, be- 
fore being driven into the linear regime. This bump in 
the values of the Piocai corresponds to a regime outside of 
the putative MW regime, when new splits in the current 
carrying channels are taking place at all length scales and 
thus there are rapidly increasing new outlets giving rise to 
the value of 2.0 for the (iiocai- New channels open, but are 
not all independent; for LxL9 arrays these newly opened 
channels will typically merge with the ground state path. 
The effective value of Piocai at 2.0 appears to be a coin- 
cidence, a malicious one for several experiments seem to 
encounter this value too. The effective exponent value at 
a given v is sensitive to the ratio of the length to widths, 
albeit in a complex fashion. 

As a consequence of finite-sizes, a crossover region over 
which the effective exponent is different from the "true 
exponent" arises. The crossover region gets larger for 
smaller system sizes. Somewhat analogous to the finite- 
size scaling exponent vt characterizing fluctuations in the 
threshold voltage, we attempt to define a finite-size ex- 
ponent vi^ which helps characterize this crossover region 
over which values of the Piocai for arrays Lx Ls deviate 
from the true exponent. From the plot in Fig. 28(a) we 
find that the best estimate is given by i^; = 1.5. Although 
the quality of the collapse is by no means satisfactory, a 
couple of trends are noticeable: there appears to be a 
a region over which the (3iocai appears to lie on a sin- 
gle curve (roughly over vL^I^^ values from 1 to 10) and 
one notices that larger systems appear to stay on the 
collapsed curve upto smaller values of vL^^^^ . It is in- 
teresting to note that it appears that vi is similar to pt, 
which if true would imply the existence of single finite- 
size length scale. Also, for small systems at large v, there 
are strong signs from the values of Piocai that the transi- 
tion towards linear behavior has begun. 

There are at least a couple of factors possibly prevent- 
ing the observance of a true scaling region. Firstly, there 
are very strong finite size effects. Systems of a sufficient 
size are required before the putative scaling behavior can 
be discerned. For example, ID channel sizes need to be 
long before the linear dependence of I on can be ob- 
served. We estimate from our analysis in Sec(Ql}that they 
should be at least longer than 1000 dots. So although 
using brute force computational power we have reduced 
significantly the statistical noise for smaller systems (e.g.. 
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FIG. 28: The collapse of the local slopes for uniform capac- 
itance systems is plotted in Fig. |2HIa) . Collapse of the local 
slopes for symmetric systems with uniform gate capacitance 
is plotted in Fig. I28f b). Note that the value of vt that gives 
good data collapse is different from the value of vt that gives 
similar collapse for systems of size L x L^''^. 



343 X 49, 216 x 36) , these systems sizes are insufScient 
to actually observe the putative scaling and we observe 
an effective exponent not in agreement with the theoret- 
ically expected scaling values. Secondly, statistical noise 
needs to be reduced significantly further for larger sys- 
tems. The reduction of statistical noise for large system 
sizes, especially at lower v, is strongly dependent on the 
cost of determining correctly the value of the current. We 
will elaborate on this further. 

Additional complications arise from the fact that there 
are large fluctuations and large timescales associated 
with channel formation. It is the complexity associated 
with both determining correctly the channel structure as 
well as the converged current value that makes the exact 



and proper simulation of electron flow in arrays such a 
difficult task. The two issues are in someways aspects 
of the same problem - the timescales required to form 
a steady state current pattern are long and broadly dis- 
tributed between samples. This phenomenon is common 
to several other dynamical systems involving collective 
transport and disorderiSSiS. The timescale required for 
current patterns to reach a steady state appears to be 
different from the timescale required for the current val- 
ues to reach steady state. For a particular sample con- 
sidered, the difference in current after the last channel 
formed was only 2%, but while investigating systems of 
the same size and at similar u^s (to within a factor of 2), 
we noticed that when channels formed, there were con- 
comitant changes in the current by more than 20%. This 
wide variation is part of the problem - for it is difficult 
to estimate how much, a well formed channel will con- 
tribute to the overall current. Any adaptive algorithm 
based upon channel formation and activity isn't easy. As 
L gets larger the problem gets more acute. A somewhat 
similar problem, is the long time scale required to form 
a channel, even if there is just a single channel involved 
in conduction. 

As a consequence of the above features, it can be diffi- 
cult to determine the current value correctly. At best, we 
can strive to minimize the probability of getting an in- 
correct current value. As with other simulation schemas, 
it soon becomes a problem of optimizing a finite amount 
of resources - the reduction of statistical noise has to be 
traded off with systematic errors. Naively one would ex- 
pect that the channels that conduct most of the current 
would form early on, and thus with simulations of suffi- 
cient duration the major current carrying channels will 
have reached a steady-state, both in terms of current car- 
ried and formation. This is not necessarily the case and 
even if it were, given the broad range of times for this to 
happen between samples of a given size, it would require 
setting all simulation runs to be sufficiently long to ac- 
commodate the longest time to steady-state. In addition 
to being difficult to estimate a priori^ it would be com- 
putationally no more efficient than using an algorithm 
that determines dynamically whether the current chan- 
nels have reached a steady-state. 

After accounting for initial transient effects, we set a 
bin size to be 10000 and compute the current in the first 
two bins, based upon which we use a convergence crite- 
ria (to be described later) to determine if the current has 
reached a steady state. If the current hasn't converged, 
the bin-sizes are doubled, i.e., the number of electrons 
that we wait for to tunnel off are doubled, after which 
similar checks to determine the steady state is done on 
the next two bins. This process disregards the history 
and values of previous bins. One of the reasons this is 
done, is because it can take very long for the steady-state 
distribution to be free of initial transients and biases. It 
is difficult to use two successive bins from the same initial 
configuration for convergence to determine in a definitive 
way whether we have reached a steady state. Our ap- 



26 



proach to correctly determining steady-state current, is 
to use two different initial conditions and to simulate un- 
til they each reach values that: (i) individually converge, 
and (ii) converge with respect to each other. (This is 
somewhat analogous to the situation for simulations of 
glassy systems where at least two different starting con- 
figurations are adopted as a measure to check against 
getting stuck in a local minima while exploring state 
space). We refer to the two starting states as the "hot" 
and "cold" configurations respectively. The classification 
of hot and cold states reflects the fact that the cold initial 
configuration has been prepared by the addition of elec- 
trons so as to have a smooth spatial gradient of electron 
potential from the emitter lead to the collector lead for 
the given value of v, while the hot initial configuration 
has a smooth spatial gradient of electron potential but 
corresponding some value greater v' than the required 
value of V. 

The rate at which the local values of current change 
can be very different for hot and cold states; it is also 
typically very different for different samples. It is possi- 
ble that a simple measure of convergence like setting an 
acceptable upper limit on the percent difference between 
the values of local current in two bins before considering 
the current to have converged, mistake the slow change 
to be an incorrect convergence. Any measure of conver- 
gence whether hybrid or for a single state should take into 
account the fluctuations in the value of the local current. 
Our method for determining convergence can be summa- 
rized as follows: After reaching threshold, we initialize 
two different states - hot and cold. We start by simulat- 
ing the hot configuration until it converges after which 
we switch to the cold configuration. We check if hot and 
cold states satisfy the convergence criteria so as to distin- 
guish it from the convergence test of two successive bins 
from the same starting convergence (single convergence). 
Every time the cold configuration is checked for single 
convergence. If the cold configuration is singly converged, 
but the hybrid convergence criteria isn't satisfied then we 
switch to the previously saved hot configuration. In gen- 
eral, we check for hybrid convergence every time a single 
convergence check is performed and toggle between hot 
and cold states every time either one of them satisfies sin- 
gle convergence but the test for hybrid convergence isn't 
satisfied. This is an attempt to keep the dynamical evo- 
lution of the hot and cold simulations somewhat in phase. 
By ensuring that the individual configurations have sep- 
arately converged - once the test of hybrid convergence is 
passed, it is fair to assume that we have determined the 
steady-state current value and pattern. For the same net 
computational resource, the hybrid convergence method 
provides higher quality data™. 

Related to the ongoing analysis of statistical versus 
systematic errors, we present some final remarks about 
simulations of larger system versus more disorder aver- 
aging: just because the ability to simulate larger systems 
may exist, does not make it necessarily an efficient use 
of computational resources. Bigger may not always be 



better - for it maybe possible to get more accurate re- 
sults by simulating larger number of samples (disorder 
realizations) of smaller system sizes than smaller number 
of samples of larger system sizes. 

A system is said to be self-averaging^^ for a variable A, 
if the error in n statistically independent measurements 
of A (AA) tends to go to zero as L ^ oo, i.e., 

AA(n, L) = {< > -< A >^)/n, n > 1 (23) 

If it goes to an L-independent value the system lacks self- 
averaging. We computed AVr(n,i) for n = 10000 sam- 
ples and we find that one-dimensional arrays are strongly 
self-averaging, i.e., AVrin, L) ~ {nL''-) , where d — 1. 
Two-dimensional arrays are weakly self-averaging: in this 
case AVrin.L) ^ (nL~^/^). Generally if a system is 
self-averaging than simulations of larger system sizes is 
not counter-productive. Recapitulate that from plots of 
Piocai in one-dimensions, we saw that linear chains of 
lengths greater than a 1000 are required to observe lin- 
ear scaling; this in turn in a way sets a lower bound on 
the sizes of two-dimensional arrays. 

B. Disordered capacitance 

We plot the I-V curves and the Piocai for DC arrays 
in Fig. I29r a'l. As can be seen in Fig. EHIb) there does 
not exist a regime where a definite single value of the 
exponent describes the scaling of I with v for all system 
sizes. Unlike UC systems, we have not simulated arrays 
of size 2744 x 196 but only upto 1000x100 - which goes 
to substantiate the dependence of the putative scaling 
exponent on system sizes. 

VI. OTHER RESULTS 

We briefiy discuss results for for one-dimensional ar- 
rays with tunneling disorder. Along with the understand- 
ing of paths at threshold from section Hvl we can use it 
to infer some basic features of the ground-state path of 
RD arrays, even though we have not explicitly simulated 
2D systems with tunneling disorder. The (iiocai for ID 
arrays with tunneling disorder is shown in Fig. 1301 The 
Piocai are qualitatively and even quantitatively similar 
to ID arrays without tunneling disorder. We note that 
the value of j3iocai approaches 1 at smaller v for larger 
system sizes. The system sizes and values of v at which 
they approach 1 are essentially similar to ID arrays with- 
out tunneling disorder. So although the slow points now 
are a consequence of a combination of tunneling resis- 
tance fluctuations and local minimums in the potential 
gradient (rate differences), the basic mechanism as out- 
lined earlier for ID arrays of overcoming slow points with 
increasing voltage remains valid and thus the linear de- 
pendence on increasing v. This combined with results 
of the transverse deviation of paths at threshold, where 
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FIG. 29: I-V curves for systems with disordered capacitance 
are plotted in Fig. I29f a). Local slopes of the I-V curves are 
plotted in Fig.l^bl. 



we found that ^ scales as La , irrespective of the type of 
disorder, indicates that a priori there is no reason to ex- 
pect that spUttings will occur any differently (pre-factors 
may change) and thus the probability is very small that 
I-V scaling on introducing tunneling disorder will be any 
different in the thermodynamic limit. 

Similar to the comparisons of C, between DC and UC 
arrays, we compute the wandering exponent for RD ar- 
rays. We find that the value of the wandering exponent 
C, as shown in Fig. l^lf a^. approaches the value of | as the 
size of systems simulated gets larger. Also, as shown in 
Fig. ISir bl. the structural properties of the ground-state 
path as measured by the relationship between the depth 
and the number of end-points is similar to that of UC. 
The probability distribution of gaps and the mean lat- 
eral length of separ ation fo r gaps of size A are plotted in 
Fig. 32(a) and Fig. 32(b) respectively. Differences with 
UC if any are not signihcant. From the comparisons be- 
tween UC and RD arrays as well as UC and DC arrays, 
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FIG. 30: The Piocai for ID arrays with offset charge disorder 
and tunneling resistance disorder. The effective exponents are 
quantitatively similar to ID arrays without tunneling disor- 
der. For larger system sizes the value of Piocai approaches 1 
at smaller u. Although the slow points are a consequence of 
a combination of tunneling resistance fluctuations and small 
voltage differences, the basic mechanism of overcoming slow 
points with increasing voltage remains and thus the linear 
dependence on increasing v. 



it can be confidently said that the main features of the 
ground-state path - meandering, structure and geometry 
- are invariant to the type of underlying disorder. 

There have been suggestions based on experiments 
Ref. that the presence of the tunneling disorder could 
lead to greater transverse fluctuations because of the in- 
troduction of additional possible bottlenecks arising due 
to the large fluctuations in the tunneling resistances. 
Based upon numerical simulations, we do not notice any 
changes from the properties of UC arrays in the trans- 
verse meandering or the structure of the paths. To vali- 
date further the claim we carried out the following numer- 
ical experiment: we first computed the path in a few UC 
arrays. Keeping everything else the same, we introduced 
disorder in the tunneling resistances in the otherwise sim- 
ilar arrays and computed the paths. For the four different 
arrays we experimented with, we did not find any signif- 
icant changes in the structure of path (although actual 
values of the current densities will be different). Fig. 1331 
shows the results for one of them. Although not conclu- 
sive, this is indicative that resistance disorder at most 
changes the current density distribution for a given sam- 
ple and that change is indistinguishable when averaged 
over many samples. It is important to mention that the 
ground state paths in Fig. 1331 and similar experiments, 
were not computed using the transfer-matrix approach, 
but was dynamically determined at = 0.0. It is possi- 
ble, however, that due to greater dynamical freedom in 
selecting current flow paths at higher values of there 
still be differences in the properties of current carrying 
paths. 
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FIG. 31: Fig. |31(a)1 plots the dependence of the standard de- 
viation of "I> with L for arrays with both offset charge disorder 
and random tunnehng resistances is plotted. Note that the in 
spite of the introduction of resistance disorder the scaling ex- 
ponent of the transverse meanderings is not diffe rent fro m the 
meandering of the first path for UC arrays. Fig.js^bj' 
a comparison of the relationship between Uep and d, 



and RD arrays, 
holds for both. 
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VII. SUMMARY AND CONCLUSIONS 

Using computer simulations we can easily control the 
presence of different disorder types and thus discern the 
individual and collective effects. In doing so, we find 
that the presence of background charge disorder is the 
dominant type of disorder, and although there are some 
minor changes for arrays with variable capacitance and 
tunneling disorder, the main scaling arguments and ex- 
ponents characterizing the arrays at Vt and in the con- 
ducting regime close to Vt remain unchanged. A study 
of the interface properties in section IIIII indicated that 
the ground-state path for two-dimensional QDA should 



FIG. 32: Plots comparing the probability distribution of gap 
sizes and the dependence of mean lateral length with gap sizes 
for UC and RD arrays. The mean lateral length scales as A"? 
for RD arrays as shown in Fig. |32(b)| whereas, Fig. |32(a)| 
compares the probability distribution of the gap sizes for UC 
and RD arrays. 



belong to the same universality class as the DPRM. By 
looking at the structure and the transverse deviations of 
the ground-state path we were able to establish the con- 
nection conclusively. We also saw in sections ^ and IIVI 
that the introduction of disordered Cs does not change 
the current-scaling exponents for one-dimensional arrays 
nor of the ground-state paths. From section [V] it ap- 
pears that the scaling exponent ( for 2D arrays does not 
depend upon the types of disorder simulated either. 

The dominance of charge disorder is probably due to 
the fact that the disorder energy scale is set by the pres- 
ence of the background charge impurities, is the crucial 
energy scale of the system. This in part is due to the 
fact that the fluctuation between the charging energy of 
dots as a consequence of the particular parameter val- 
ues we choose (C™°^ = 2.0), is less than the fluctuation 
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FIG. 33: Comparing the ground state path at z/ = 0.0. 
Fig. |23b) is an array with exactly similar charge disorder as 
the array in Fig. lSSf a). but with tunneling resistance disorder 
included. Although the current densities at various locations 
are different, the overall structure is similar. 



in electrostatic energy due to offset charges being cho- 
sen randomly between [0,1[. The presence of tunneling 
disorder does not change the energetics of the arrays, 
i.e., Vt and fluctuations in Vr- It is important to re- 
mark, however, that if the offset-charge impurities were 
disregarded and only a non-uniform Cs considered, the 
arrays would still exhibit a threshold voltage, separating 
the insulating and conducting phases and most proper- 
ties would still be similar to the situation where there 
was only offset-charge impurities. If the energy fluctu- 
ations due to of non-uniform Cs, was greater than the 
background charges the claim would be that the domi- 
nant form of disorder was the Cs although the properties 
of the array would be essentially insensitive to which was 



the dominant disorder. 

From our discussion in section 0, we make the im- 
portant conclusion that it is most likely that the MW 
hypothesis is correct and valid for disordered QDA irre- 
spective of the actual relative strengths of the disorder. 
It is important to appreciate that one needs to get suffi- 
ciently close to threshold to observe the scaling and that 
too only for large systems. 

From our experience, a naive approach to determining 
the steady state current consistently underestimated the 
current values, which tends to get more acute at lower 
values of v. As a consequence, a higher putative value of 
(3 would be observed. It also follows that there is a need 
for careful simulations, for with slightly less diligence it 
would have been tempting to predict an exponent range 
of 2.0-2.25. This is intertwined with the issue of high 
computational cost, which arises from a combination of 
the need to compute the converged current accurately for 
a single sample and the need to simulate large systems as 
a consequence of strong finite-size effects. It is not easy 
to formulate an elegant algorithmic solution to this prob- 
lem. Although, parallelization is a well defined and often 
used approach to reduce time-to-solution of a problem, 
our problem does not appear to be a suitable candidate, 
for as mentioned, one of the primary bottlenecks in our 
simulations is the extremely long times required to reach 
a steady-state configuration. It is physically-meaningless 
to run a simulation at a time t2 without state informa- 
tion at time ti, where t2 is a time later than time ti; thus 
there is a fundamental limitation on speed-up that can 
be achieved via parallelization. But parallelization pos- 
sibly along the lines of Ref. 'iS'l may be a possible route 
forward. 

We have focused on QDA in the extreme limit 
where the screening-length is less than the spacing be- 
tween dots. The opposite regime of essentially infinite 
screening-length has been well studied, both numeri- 
cally and theoretically for non-disordered arrays^ and 
recently for arrays with a random background potentiaU^ 
- although neither of these studies, nor others that we 
are aware of, use the statistical physics approach that we 
have used. Surprisingly, there has been little activity in 
the regime representing the middle ground, viz., a screen- 
ing length from a few upto a dozen dot spacings. With 
a screening-length more than a single dot spacing, the 
on-site interaction model that we have used in this work 
is not valid and computational approaches will require 
fundamental reworking. Ironically this regime is impor- 
tant (and interesting), as most nanoparticle arrays as a 
consequence of the absence of an underlying gate most 
probably have an electrostatic screening-length of a few 
dot spacings. 

In summary, we have investigated the effect of disor- 
der on the transport of electrons in arrays of mesoscopic 
sized metallic islands, at, below and above a critical volt- 
age Vt. In contrast to experiments, using computer sim- 
ulations we can easily control the effects of different dis- 
order types. We find that the presence of background 
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charge disorder is the dominant type of disorder and al- 
though there are some minor changes with the addition 
of variable capacitance and tunneling disorder, the main 
scaling arguments and exponents characterizing the ar- 
rays at threshold and in the conducting regime remain 
unchanged. Our numerical results indicate a value for 
the exponent [3 to be in the range 1.5-2.0. 
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It can be argued that if we are simulating systems of sizes 

2 

L X La , then the wandering exponent can be no larger 
than |. To establish that C, is not system width limited, we 
simulated systems^ with multiple basins, i.e., with widths 

2 

greater than L 3 . We find from a plot of the effective local 
exponent values for C,, that the value of C, is still consistent 
with I 



